FUSION OF THRESHOLDING RULES DURING WAVELET-BASED NOISY IMAGE COMPRESSION

The new method for combining semisoft thresholding rules during wavelet-based data compression of images with multiplicative noise is suggested. The method chooses the best thresholding rule and the threshold value using the proposed criteria which provide the best nonlinear approximations and take into consideration errors of quantization. The results of computer modeling have shown that the suggested method provides relatively good image quality after restoration in the sense of some criteria such as PSNR, SSIM, etc.


Introduction
Images which are registered from different sources and which have to be transferred or archived can be distorted by specific noises having a multiplicative character.In many cases, video sources and data links add their own noises during image formation and transferring.For example, we can discuss about synthetic aperture radar (SAR) images with multiplicative noise known as speckle, infrared devices with fixed pattern noise (FPN) like unstable photo element's voltage sensitivity, etc.There are a lot of methods and algorithms invented in last decades trying to effectively filter and/or compress noisy images.In the most cases, the tasks of filtering and compression are decided separately.Therefore, filtering and compression algorithms are not harmonized, that leads to new distortions and artifacts onto the images restored after compression.Moreover, being theoretically focused on an additive Gaussian noise model, any filtering algorithm inevitably leads to unsatisfactory results.
The good idea is to combine filtering and compression within one and same procedures of noisy image processing.The main decisions which we have chosen for this overview are relatively old and can be grouped into three sets.The first set of methods presented by works [1,2] was developed using the Rissanen's principle of minimal description length (MDL) [3].Here, the final compression ratio is determined by the achieved maximal quality of the decompressed image.The second set of methods presented by works [4,5] exploits the idea of so-called optimal operation point (OOP) which is the compression ratio or bit rate (BR) measured in bit per pixel (bpp), when the maximal peak-to-signal noise ratio (PSNR) is achieved.Existence of OOP can be explained by existence of so-called "dead zone", an interval of quantization near zero, size of which is changed in dependency on the given bit rate for any codec.Again, OOP provides the desired quality of compressed images but does not guarantee the necessary bit rate for transmission through channels with fixed bandwidths.The third set of methods presented by our works [6,7] and other works is based on forming an image of high quality under the given bit rate (quota of bits).
The most of mentioned works uses wavelet transformation which allows consolidating filtering and compression by thresholding techniques.Usually, so-called fast discrete wavelet transform (FWT) [8,9] is the base of many suggested methods and algorithms.This wavelet transform in the form of dyadic decomposition has been applied to image processing by S.Mallat [8].The image is decomposed into several high-frequency pseudo images (subbands) containing wavelet coefficients of details in horizontal, vertical and diagonal orientations with increasing scale.Different thresholding rules are the bases of wavelet filtering methods [10,11].There are some popular shrinkage policies: semisoft, non-negative and n-degree garrote and hyperbole thresholding [10].Nevertheless, semisoft thresholding where hard, soft and Vidacovic thresholding rules are related is the most popular in practice because the one parameter (the threshold value) is needed to calculate only.Otherwise, non-negative garrote rule [10], for instance, requires two parameters to be estimated that gives an additional optimization problem.Concentrating on semisoft thresholding rules in this paper, we can note that the effectiveness of filtering for any thresholding rule depends on different image features such as texture, type and intensity of noise, spatial and radiometric resolutions, etc.In other words, each of thresholding rules can win in the contest for the best quality of the filtered image.Moreover, the situation becomes unpredictable when the noised image is undergone by wavelet-based compression that gives errors of so-called non-linear approximation and quantization [12].Therefore, this paper comprises an attempt to develop the theoretical base to fuse different semisoft thresholding rules during data compression of noisy images that provides the fully automatic scheme of choosing the best thresholding rule and the corresponding threshold value.
This paper consists of five sections.After this introduction (Section 1), Section 2 describes the problem definition.Section 3 contains the obtained criteria for different thresholding rules and discusses what to expect from the fusion of the thresholding rules.The results of computer modeling and comparison are represented in Section 4. Conclusion remarks are reflected in Section 5.

Problem definition
Let us assume that the compression is subjected to the observed signal Y, which is the projection of the unknown original image X distorted by multiplicative noise Z with unity mean and variance
(1) The noise Z may have different probability density function (pdf), for instance, an exponential pdf for speckle in SAR images; normally distributed for FPN in infrared sensors, etc. Multiscale analysis provides a decomposition of a noisy signal (1) in the form of FWT for a given number of levels Q.Because the wavelet transform is the result of successive convolutions, we define an operator .Then, the wavelet decomposition of a noisy signal with a fixed basis may be represented as follows: , where Therefore, in wavelet domain, the multiplicative model ( 1) becomes the additive one (2).
Any wavelet coefficient can be represented as ( j I is the quantity of wavelet coefficients at each level j, Q j ,..., 1 = , and the noise in wavelet coefficients can be shown as: , and here, we consider conditionally that i e is Gaussian noise with zero mean, and .Therefore, we complicated the model for noise which is used for computer modeling and reflects heterogeneous features of an image in FWT subbands.There are some works where different hypotheses about noisy wavelet coefficients distribution are suggested.The generalized Gaussian distribution has been very well adopted in imagery on wavelets (see, e.g., [13]).Here, we consider that the high-pass filtering leads to an approximately symmetrical form of the distribution for wavelet coefficients.Hence, the properties of multiplicative noise in spatial domain are changed in such manner that we can use the additive white noise model in wavelet domain.This is our lonely presumption which has no the strong mathematical background.
Thus, the problem is to perform data compression of a noisy signal (1) by encoding the wavelet coefficients (3) in such way that we could improve and estimate a quality of the processed noised image.It can be done by use of mean square error (MSE) criterion in the form of the Euclidean norm which can be written as follows where X W ˆ is the wavelet coefficients of the restored (after compression) image.
It is because the original signal is distorted, then the coefficients X W ˆ should be considered as estimators of the wavelet coefficients for the original signal X after processing of the noisy wavelet coefficients Y W .It is shown in the papers [8,9] that the error of restoration (4), which is computed in the wavelet domain, is equivalent to the averaged square of the norm of the error in the spatial domain: Let us suppose that the noisy signal (1) is represented by I discrete samples (pixels); in the case of FWT, it gives the same quantity of the wavelet coefficients I.Then, after compression, we have M so-called significant wavelet coefficients only ( I M < ) which are undergone by quantization before statistical encoding.If we ignore the possible small losses of the signal energy during encoding, then the estimators of the wavelet coefficients can be represented as follows: where the brackets ⋅ denote the operation of uniform quantization (uq).To calculate the total error caused by the approximation and quantization, we first consider the sum of squares of the errors in estimating the wavelet coefficients: For the sake of notation simplicity, we enter the next terms: is the variance calculated for all of wavelet coefficients; is the variance of the thresholded (significant) wavelet coefficients, is the variance of wavelet-coefficients for noise with zero mean; is the quantization variance for M significant wavelet-coefficients.
Then, the formula for the MSE criterion can be written as follows It follows from equality 8 that minimizing the value is equivalent to maximizing the value ( ) . In order to obtain the best MSE, it is necessary to provide the maximal variance of the significant wavelet coefficients (because it gives the good signal approximation [8]) and the minimal variance of quantization errors.These requirements are contradictory, because the larger M of the significant wavelet coefficients increases the total quantization error simultaneously.On the other hand, we can control both these components in equality 8 by different methods including thresholding techniques.Therefore, our problem can be described in the terms of the approximation theory where any thresholding rule can be considered as a procedure of non-linear approximation of the noised signal (image) when the maximal value ( ) thresholding rules, we have a nice opportunity to choose the best rule because we obtain different signal approximations.Therefore, the problem converges to the search of the generalized method suitable to estimate the effectiveness of any thresholding rule.Our solution placed below has been derived for so-called semisoft thresholding rules which are described in Section 3.For computer modeling and comparison, the quality of the restored after compression image is estimated by several criteria such as peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM).

The fusion of thresholding rules
The main idea of our method is to iteratively refine a non-linear approximation of the noised image using so-called coherent structures which are the components of the image and which are strongly correlated with a chosen wavelet basis [15].Having a variance for any orthonormal wavelet basis.The parameter I ρ can be considered as the supreme estimator of the correlation coefficient for correlation between Gaussian noise and any given wavelet basis.It follows from equality 9 that this estimator does not depend on noise variance.It has been experimentally proven by our experiments [6] that the theoretical value of I ρ is also the highest asymptote of the correlation coefficients for correlation between Daubechies wavelets, symmlets, some biorthogonal wavelets, on the one side, and, on the other side, noises having different pdfs and their combinations, e.g., log-normal, exponential, Gamma-distributed, etc.
It is known from the theory of wavelets [8,14] that the nonlinear approximation of any signal decomposed by wavelet transform is referred as inverse wavelet transform W -1 for the first (significant) M wavelet coefficients after their sorting ˆ. (10) where The residual image M Y can be computed as the difference between the input image and the pseudo image obtained by means of the inverse wavelet transform for the coherent structures The residual image M Y can not be recognized as a noise if the following inequality 2 ) ( (13) is true.Not taking into consideration the quantization errors, the estimator X ˆ of the original image is the sum of M coherent structures (10).
Hard thresholding.The equality 13 shows that pursuit of the coherent structures looks like hard thresholding (see Figure 1, a) of the wavelet coefficients using the threshold [8,10,14]

IME T &
Therefore, the threshold value τ is automatically set by any change in the quantity M of the significant wavelet coefficients.Taking into consideration equality 8 and the following conclusion from Section 2, we can derive the method of the noised image processing based on hard thresholding: ( ) Let the wavelet-coefficients of the noised image be arranged according to equality (11).Then, the proposed algorithm using the matching pursuit can be represented as follows [6].
Compute the cumulative sum of wavelet-coefficient squares and, according to (14) calculate the threshold value for the given correlation coefficient Otherwise, go to the step 2. It is possible to use different rules to finish calculations when the maximum of the cost function ( 15) is achieved, for example , where δ is the reasonable error of calculations showing the accuracy of the algorithm; In practice, the processing scheme in the form ( 17) is applied to the details only; while the wavelet coefficients of the approximation remain unchanged.Therefore, we need to analyze the three components.
where t is the number of the first-ordered wavelet coefficients corresponding the approximation.The third (last) term in equality 18 represents the error which is inevitably appeared due to the thresholding of the wavelet coefficients.We obtain for the first two summands in equality 18:  ) and (Δ is an interval of quantization).Finally, the equality 8 is modified for the case of soft thresholding: Hence, we need to also modify the method of processing based on equality 15 and equality 16: The algorithm corresponding to equalities 21 and 22 contains the same steps with some modifications as it is for the case of hard thresholding.The estimators of the processed wavelet coefficients can be calculated as follows: , .
Vidacovic thresholding.The thresholding function derived by Vidacovic [10] comprises some "medium" variant of semisoft thresholding as it can be seen from Figure 1, c: We obtain after some algebra: The corresponding algorithm for noisy image processing can be also described by two equations:

IME T &
The estimators of the processed wavelet coefficients can be found as follows:   Fusion of thresholding rules.In order to fuse different thresholding rules, we need to perform some formalization.Suppose that there is a set of q thresholding algorithms Therefore, equality 29 determines the best thresholding algorithm * i A (or/and the corresponding thresholding rule).
The obtained wavelet coefficients are sent to any wavelet coder where they are statistically encoded and transferred.The restored image can be considered as the sum of the coherent structures (in the form of equality 10) which are restored by the corresponding decoder on the side of recipient or from the database (in the case of archiving).
Our expectations.In the case of data compression of noisy images, we have too many factors influencing on the final result.Therefore, we need to determine the strategy of using the suggested method which can answer, by our opinion, on the following two important questions: 1) What, in average, have we to expect from the fusion of thresholding rules?2) How "to tune" the suggested method for practical applications, when, for example, bandwidths of data links are limited?In order to obtain an answer for the first question, we should remind that there were a lot of theoretical works devoting to wavelet shrinkage.In last decades, Mallat, Coifman, DeVore, Donoho, Johnston and their collaborators explored statistical optimality and risks of wavelet shrinkage for cases of signal processing, approximation, and statistics (see, e.g., [8,10,15,16]).
These works help to find out the high estimators of risks for any thresholding rules; however, it does not allow us to compare or choose the best thresholding rule in the given cases.

IME T &
On the other side, there were a lot of works containing the results of practical applying of the semisoft thresholdiong rules to different fields of imagery.Unfortunately, in those works, comparison is limited by types of images and noises (see, e.g.[17,18]).If we take this way for proving effectiveness of the suggested method, then we can put in the difficult situation of endless comparing by pairs, for example, hard and soft thresholding under different bitrates, etc.
Nevertheless, here, we try to give our (maybe, "too engineering") answer for the first question; our thoughts are not controversial to the main theoretical results.
When we consider a noisy image, then the behavior of the ordered wavelet coefficients (equality 11) is changed; we observe a difference between two curves for small wavelet coefficients (details of the FWT) (see Figure 2).Differences depend on intensity of noise.It means that different noise variances lead to different threshold values; however, it does not mean that any thresholding rule will be a winner in the given case.Figure 3 (in the form of sketch for the case of a strong noise) shows some different situations which can be encountered during modification of wavelet coefficients in accordance with thresholding rules ( 14), (17), and (24).These situations have the following explanations.
If the input image ( 1) is relatively clear (noise is too small) then we observe that the quantity of undistorted wavelet coefficients with small amplitudes (details) is relatively large; hence, we can choose the small threshold value M Y w to set these coefficients to zero for compression.It is obviously, that the hard threshodilng rule has to win because no modification is done in the residual wavelet coefficients.In the case of noisy image, the situation is not so obvious because the dynamic rate of wavelet coefficients related to details is increased due to noise; hence, the threshold value should be increased.If we register the case of a strong noise, then, we consider, the soft trhesholding rule should be a winner because the residual distorted wavelet coefficients are changed (smoothed).It can be noted from Figure 3, where "soft" modification sets the part of the significant wavelet coefficients (the residual details) to be close to "an ideal", high frequency noise-free coefficients.Also, as it seen from Figure 3, in general, the Vidacovic thresholding rule takes an intermediate position between soft and hard threshoding rules.It means that the Vidacovic thresholding rule modifies wavelet coefficients related to middle frequencies of the signal.Hence, we cannot predict all possible results.
Therefore, we should go to the second question because error of quantization depends also on available bit rate for data links or the given quota of bits for image archiving.The problem is that the interval of uniform quantization depends on both the given bit rate and the dynamic rate of the residual wavelet coefficients: where C R is the given budget of bits; 1 Y w is the first coefficient from the ordered wavelet coefficients having the maximal value.In the cases of low bitrates ( C R <I), there are possibilities to tune a full scheme of noisy image compression if to determine so-called optimal operation point (OOP) [4,5] of a coder.It has also been experimentally proved by the author for a case of multiplicative noise that any coder has such bit rate (i.e.OOP) where the maximal peak signal-to-noise ratio (PSNR) is achieved.Therefore, our experiments have to show that the suggested method chooses different thresholding rules and OOPs in dependency on noise variance and image textures.

Results of modeling
For our experiments, we formed the library of test images containing 25 reference images which have been taken from the library described in [20].The initial formats all of images were converted into the gray scale format with 8 bpp.Moreover, then, all of the converted images were divided into two groups with homogeneous and heterogeneous textures.The type of texture was determined by use of so-called variation coefficients which can serve as indicators of homogeneity within the given window [17].Image is considered to be homogeneous if the quantity of points (pixels) from the set of homogeneous points is more than the reasonable threshold.We set this threshold value to 60 %.We used the FWT of three levels (Q=3) based on the wavelet CDF 9.7 (Cohen-Daubechies-Feauveau).To simulate a multiplicative noise, the random number generators were applied to model speckle-noise with exponential pdf and unity mean.The variance of noise was a subject to change during experiments.In order to get a completely full processing of noised images, we used SPIHT (Set Partitioning In Hierarchical Trees) coder software [19] which has been modified to regulate the dead zone in dependency of the was changed within low bitrates (up to 1 bpp) with the step 0.01 bpp in order to find out the OOP.
The results of computer modeling are represented in Table 1.The table contains results in the form of different characteristics such as PSNR, SSIM [21] and OOP.The results obtained by the suggested method (i.e., the chosen thresholding rule) are marked.We compared the obtained results with "not-chosen" thresholding rules and SURE and Oracle thresholding rules under the same OOP.The SURE thresholding algorithm is described in [8]."An ideal denoising method" was represented by the principle of "Oracle" when the threshold value is calculated if the original (non-noised) wavelet coefficients are known [8,16].
Some distorted and restored after compression images with the obtained optimal thresholding rule and bit rate are shown on Figures 4-11.One can see from Table 1 and Figures 4-11 that the suggested method provides better image enhancement both numerically and visually.
In spite that there are other thresholding rules, we do not place here results of comparison with these rules because it is difficult to find out any statistical lows during pair comparing.Moreover, we understand that our method cannot be a winner in many cases if to compare to powerful and very complex methods as [1,2,6,7,13,14,22].

Conclusion
The problem of noisy image compression is very complex; and, certainly, the suggested method does not pretended to decide related tasks in full.The main aim of this paper is to demonstrate that there are possibilities to build more flexible schemes of noisy image processing.The suggested method can be expanded by means including other thresholding rules, e.g.non-negative and n-degree garrote.However, it will require using a multi-criteria problem definition.
Because the suggested method exploits the idea of pursuit the coherent structures, than it leads to relatively high computational and timing expenditures.We estimate the complexity of our method as three thresholding rules and choosing the quantity M of the significant wavelet coefficients is similar to search of the best path during wavelet package processing.In spite, we consider that there are some possibilities to decrease these expenditures due to using socalled pipe-line and unrolling loop methods of programming to seek the maxima of cost functions (15), (21), and (26) in parallel.It allows us to promote this method for onboard mage processing using FPGA implementations.

3 .
Compute the new (corrected) value M of the significant wavelet coefficients from the inequality (16).4. Check: is the maximal value of the cost function (15) achieved? 5.If the cost function(15) has the maximal value then calculations are finished.

M
is the estimator of the quantity of the significant wavelet coefficients obtained at the iteration with number n, n=2, 3, ... Because the algorithm presented by equalities 15 and 16 is based on the hard thresholding rule, then the estimators of the processed wavelet coefficients are automatically obtained as Soft thresholding.During soft thresholding (see Figure1,b), the wavelet coefficients are also set to zero if their values are below the threshold, but simultaneously, there is a reduction of the details (the wavelet coefficients of high-frequency subbands of the FWT) to the threshold value showing how many wavelet coefficients of the approximation are in the total quantity of the wavelet coefficients.We assumed during our derivation i

,
which can search for coherent structures in the chosen wavelet basis.Let the results of executing the corresponding algorithms be the values of the objective functions select the best coherent structures, we must use the criterion:

4 ( 3
log ) O II operations taking into consideration that we have DOI 2 )

Table 1
The results of comparison.