The use of non-linear adaptive filters in noise removal was discussed in 4.6. In this section, we will show how the discrete wavelet transform can be used to remove noise from a corrupted image using a process known as •wavelet denoising. This technique, initially proposed by Donoho,24’25’26 relies on the idea of thresholding in the wavelet domain and remains very much in vogue and the subject of active research. The accolades attributed to wavelet denoising are certainly impressive – advocates claim that it "offers all that we might desire of a technique, from optimality to generality."27 While there may be some amount of hyperbole as to it’s efficacy, Donoho did show in his seminal work that the technique attains near-optimal noise reduction properties, in the minimax sense, for a large class of signals – although it should be stressed that much of this work assumed the corrupted signal suffered from additive Gaussian noise. One of the driving factors behind use of the technique is its inherent simplicity. For these reasons, after image compression, wavelet denoising represents perhaps the second most common class of wavelet-based algorithms with regards to image processing. In this section, the theory behind wavelet denoising is introduced, a MATLAB prototype is discussed, and the use of the Wavelet Toolbox for image denoising is also briefly touched upon. Following these preliminaries, we take a quick detour into explaining how to implement an efficient 2D IDWT algorithm in C. Wavelet denoising requires the inverse transform and as we shall soon see, both the C62x/C67x and C64x versions of IMGLIB do not provide inverse wavelet transform functions. Once we have C functions for the 2D IDWT, we can proceed to implement wavelet denoising on the C6416.
Consider an image f(x,y) corrupted with additive noise n(x,y). Then the observed image g(x,y) is
We wish to recover from the noisy image g an approximation f to the original f, which minimizes the mean squared error between the approximation and f. Let W[-] and W"’[-] denote the 2D DWT and 2D IDWT, respectively. Now let T[-,X] be a wavelet thresholding operator with threshold X. Given these definitions, a wavelet based denoising procedure is summarized in Algorithm 6-2.
Algorithm 6-2: Wavelet denoising
Although some researchers have reported slightly better results using the overcomplete wavelet transform we used in the preceding section for multi-scale edge detection,28’29 the typical means of implementing wavelet denoising of images is to use the classical 2D DWT, as described in 6.1. The algorithme a trous has comparatively larger storage requirements due to the use of the undecimated transform, and concomitantly more comparisons needed for the thresholding operator. These two facets combine for a substantial increase in computations that limit the practicality of utilizing the undecimated wavelet transform in image denoising algorithms.
Thresholding in the wavelet domain works to remove noise due to the "energy compaction" property of the wavelet transform. This property is an integral part of what makes the DWT such a useful tool in signal processing, and while we apply it here in the context of image denoising it is equally vital in image compression. What this property means is that as we proceed with the wavelet transform, signal energy concentrates in the LL band, while the noise energy dissipates throughout the other detail subbands. Hence we use the DWT to separate the actual signal content from the noise, where the low-frequency "base" of the image is contained in the LL wavelet coefficients. If the coefficients in the detail subbands are significant (large in magnitude), they are attributed to the original signal.
![Hard threshold versus soft threshold using input data normalized to the interval [-1,1]. (a) Hard threshold transfer function, with(b) Soft threshold transfer function, with](http://what-when-how.com/wp-content/uploads/2011/09/tmp17F186_thumb.png "Hard threshold versus soft threshold using input data normalized to the interval [-1,1]. (a) Hard threshold transfer function, with(b) Soft threshold transfer function, with")
Figure 6-16. Hard threshold versus soft threshold using input data normalized to the interval [-1,1]. (a) Hard threshold transfer function, with
(b) Soft threshold transfer function, with
We exploited this fact in the multiscale edge detection algorithms of the previous section, as we saw that in a multi-level wavelet decomposition, the number of wavelet coefficients with significant magnitude is small, and the significant coefficients correspond to high-frequency detail and edges. Insignificant detail coefficients are assumed to be noise and by discarding these coefficients corresponding to n(x,y) we are hopefully able to reconstruct j[x,y) without loss of detail.
This procedure works so long as the noise power is much smaller than the signal power, the signal is assumed to have low frequency components, and n(x,y) is white1. Although the above description may make it appear that it is a form of smoothing, wavelet denoising is non-linear and not a smoothing operator. A smoothing operator destroys high-frequency content at the expense of low-frequency content, while wavelet denoising "attempts to remove whatever noise is present and retain whatever signal is present regardless of the frequency content of the signal."27
Algorithm 6-2 leaves out two important facets to the denoising algorithm: the details of the thresholding operator T and selection of the threshold X. Generally speaking, T[-,A,] comes in two flavors, a so-called hard threshold or soft threshold, examples of which are plotted in Figure 616. Hard thresholding of wavelet coefficients is a form of thresholding similar in form to that described in 5.2.1 and implemented in various functions within IMGLIB and the Intel IPP library. This operator retains all coefficients greater and than A, and sets the remaining coefficients to zero:
Inspection of Figure 6-16a shows an abrupt discontinuity in the hard threshold transfer function at |x| = X. This discontinuity tends to manifest itself as annoying artifacts in the denoised image, and motivated researchers to propose the soft threshold scheme, otherwise known as "shrinkage". This transfer function is more continuous and aims to shrinks those coefficients above X in absolute value. The soft thresholding operator sets to zero all wavelet coefficients x if |x| < X, and otherwise:
In general, there are three different methodologies for performing the threshold operation on the DWT coefficients. Most of these do not apply any thresholding on the LL (approximation) subimage:
1. Universal Threshold: a single threshold is applied globally to all detail coefficients, at all scales. Donoho and Johnstone used this methodology in [25] as part of their VisuShrink wavelet shrinkage denoising algorithm. This threshold has been shown to be a good choice for use with soft-thresholding on images corrupted with additive Gaussian noise, although it does tend to produce over-smoothed images.
2. Level-Dependent Threshold: The SureShrink16 and BayesShrink30 methods attempt to improve upon a single threshold X by making it subband-adaptive. A different X is computed for each detail subband in each wavelet scale.
3. Spatially Adaptive Threshold: Akin to the adaptive filtering schemes of 4.6, these operators build upon (2) by tailoring X based on the local characteristics of the subband. By their very nature, they usually scale-dependent.
Determination of the actual value(s) used for X is critical to the performance of wavelet denoising, and is not an easy task. Here we find ourselves in the same boat as in next topic, where it is highly desirable to have a data-driven means of calculating X. Many schemes abound, most of which compute X based on an estimate of the signal energy and noise variance a2 (see 4.6.1 for an explanation of this parameter). Donoho and Johnstone proposed in [25] as part of VisuShrink the universal threshold
where M is the sample size, or total number of image pixels in our case. VisuShrink uses a robust estimation of a based on the median absolute value of the DWT coefficients:
where Wc denotes the wavelet coefficients. While we use this estimate of the noise power in the MATLAB prototypes that follow, it is difficult to make a case for it in an embedded DSP wavelet denoising implementation. DSPs may excel at certain numerical computations; however they begin to lose some of their luster when utilized in more logic-based operations. Sorting is one such operation, and in order to use the median estimator we must sort a potentially rather large set of numbers so as to calculate the median wavelet coefficient value. The partial sorting optimization described in 4.5.4.3 does not apply here because we are not sorting a fixed, small neighborhood of values. There are parallel sorting algorithms that are designed to sort large arrays of data in chunks. Some variant of this procedure would be needed to efficiently sort the DWT coefficients. The best means would be to split the data into blocks small enough to fit in internal RAM, sort each of these blocks, and then merge the results into a single sorted array.
A popular level-dependent variation on the universal threshold theme is:
again where M is the total number of signal samples and j is the wavelet scale level31. The algorithm we implement on the DSP is a subband-adaptive, hard thresholding denoising algorithm that uses the standard deviation of one of the detail subbands (HH) as an estimate of the noise power. This procedure is summarized in Algorithm 6-3.
Algorithm 6-3: Level dependent hard threshold wavelet denoising INPUT: Noisy image I, # decomposition levels N, multiplier s OUTPUT: Denoised image K
Wavelet Denoising in MATLAB
The Image Processing Toolbox function imnoise has been introduced in previous topics as a convenient means of generating test data for image processing simulations. If this toolbox is available, the following commands can be used to generate pixel data suitable for image enhancement via wavelet denoising:
This MATLAB code snippet contaminates the image I with zero-mean Gaussian white noise of a specific noise variance, determined by the variable sigma. If the Image Processing Toolbox is not available, the M-file add_zero_mean_gwn.m, located in the Chap6 directory on the CD-ROM, can be used as a replacement. This M-file uses the built-in MATLAB function randn, which returns zero mean Gaussian white noise with unit variance. The input pixel data is scaled appropriately, the output from randn multiplied by the square root of cr, and finally the two are added together before the returned data is casted back to the original MATLAB class type. This function, shown in Listing 6-13, was used to generate all of the test data for the wavelet denoising simulations presented in this section.
Listing 6-13: MATLAB function for generating zero-mean Gaussian white noise of a specific magnitude.
Two MATLAB wavelet denoising prototypes are located in the Chap6 directory. Both functions use the multi-level D4 system (see 6.2) for the DWT and draw on some of the functions introduced and developed in 6.1.1. The MATLAB 2D DWT code in 6.1.1 utilized a function haar basis (see Listing 6-1) that returned the polyphase matrix H, constructed in such a manner that convolution followed by downsampling by 2 was achieved via matrix multiplication of H and the input signal. The D4 wavelet transform is defined by the following four tap low-pass filter4:
Hence a similar periodized polyphase matrix H is needed that contains the circular shifted D4 filter coefficients, in order to utilize a linear algebra based implementation of the DWT. This matrix will be of the form explained in 6.1.1, and the MATLAB function daub4_basis (Listing 614) creates this matrix, which is then used to implement the forward and inverse D4 transform.
Listing 6-14: MATLAB function for the creation of the D4 forward transform matrix.
Figure 6-17. Wavelet denoising with the D4 wavelet on the "Barbara" image corrupted with additive Gaussian white noise. Images (c) and (d) were generated via denoise_daub4_universal, using a three-level wavelet transform and threshold multiplier of 0.5. (a) Original image, (b) Corrupted image with a noise variance of .005, created using the add_zero_mean_gwn function from Listing 6-13. (c) Hard thresholding result, with a threshold derived using the universal threshold, (d) Soft thresholding result, using the same threshold value as (c).
Both MATLAB prototypes implement hard and soft thresholding of DWT coefficients. In denoise_daub4_universal, a global threshold similar in concept to VisuShrink is used to perform wavelet denoising. The value of X is calculated from the finest HH subband, using the median estimator of the noise variance. Clients of this function may pass in a scalar multiplier (which defaults to 1) in order to adjust this estimate, if so desired. Furthermore, the function accepts the number of wavelet decomposition levels and a threshold type string (either ‘H’ for hard thresholding or ‘S’ for soft thresholding). This function is given in Listing 6-15 and some results using the "Barbara" image are shown in Figure 6-17.
Listing 6-15: MATLAB M-file for multi-level wavelet denoising using a global threshold and the D4 wavelet.
The denoise_daub4_subband_adaptive .m M-file is a subband-adaptive wavelet denoising function that implements the scheme mapped out in Algorithm 6-3. This function uses the standard deviation of each HH subband to derive X for the current level. The built-in MATLAB function std2, which returns the standard deviation of a matrix of values, is used as a convenient means of calculating the threshold value. Results for the "Elaine" image are shown in Figure 6-18, and relevant portions of the M-file are given in Listing 6-16. It should be noted that both Listings 6-16 and 6-17 leave plenty of room for optimization. As we shall soon see when we take a closer look at an efficient C implementation of the IDWT, sub-optimal wavelet implementations perform superfluous multiplications by zero, and our MATLAB linear algebra based approach suffers from the same problem. The DWT and IDWT are actually examples of multi-rate signal processing and there are functions available in the Signal Processing Toolbox that can be used to efficiently pass signals through QMF banks. The upfirdn32 function is just what the doctor ordered when it comes to implementing the IDWT in an optimized fashion. This function upsamples, FIR filters, and then downsamples a signal, precisely the operations we need to perform the IDWT. However the goal of these M-files is not to illustrate an optimized MATLAB-centric means of performing the DWT and IDWT, rather they are meant to provide a fairly easy-to-read reference implementation of Algorithm 6-3. As we have previously seen, a straightforward reference implementation is worth its weight in gold during development and deployment to an embedded platform.
Listing 6-16: Portions of a MATLAB subband-adaptive wavelet denoising implementation. See Listing 6-15 for the contents of the hard threshold and soft_threshold functions.
The Wavelet Toolbox includes a gamut of useful functions for wavelet denoising of signals and images. These range from catch-all functions like wdencmp that perform everything pertinent to denoising (controlled via a list of string input arguments) to lower-level functions such as wnoisest (noise estimation) and wthresh (soft or hard thresholding). By piecing together various Wavelet Toolbox functions it is possible to very easily prototype any number of denoising schemes. In addition, the Wavelet Toolbox includes a GUI program that puts a user-friendly face onto the myriad number of wavelet denoising options available through the toolbox. Interested readers should consult [6] for further in-depth information regarding these functions and the GUI application.
