## Denoising: wavelet thresholding

We have seen a naïve way to perform denoising on images by a simple convolution with Gaussian kernels. In this post I would like to compare that simple method with a more elaborate one: *wavelet thresholding*. The basic philosophy is based upon the local knowledge that wavelet coefficients offer us: Intuitively, small wavelet coefficients are dominated by noise, while wavelet coefficients with a large absolute value carry more signal information than noise. Being that the case, it is reasonable to obtain a fast denoising of a given image if we perform two basic operations:

- Eliminate in the wavelet representation those elements with small coefficients, and
- Decrease the impact of elements with large coefficients.

In mathematical terms, all we are doing is *thresholding* the absolute value of wavelet coefficients by an appropriate function. For example, if we chose to eliminate all coefficients with absolute value less than a given threshold , and keep the rest of the coefficients untouched, we end up with the thresholding function given by:

We refer to this function as a *hard thresholding*. But we might want to decrease the impact of the elements with large coefficients in absolute value—this is what we refer as *soft thresholding*. A common way to perform this thresholding is by simple imposition of continuity of the function; e.g.:

Other thresholding functions are of course possible, and the reader is encouraged to come up with different examples, and experiment with them.

The wavelet module `PyWavelet` is easily embedded in `sage`, and allows us to perform very fast and elegant denoising codes. The following session shows an example: We load the image of Lena as an array of size containing values between `0` and `255`. We contaminate the image with white noise and perform denoising with this technique. Since we require all nine levels of the Haar wavelet coefficients of the image.

from numpy import * import scipy import pywt image = scipy.misc.lena().astype(float32) noiseSigma = 16.0 image += random.normal(0, noiseSigma, size=image.shape) wavelet = pywt.Wavelet('haar') levels = Integer( floor( log2(image.shape[0]) ) ) WaveletCoeffs = pywt.wavedec2( image, wavelet, level=levels)

The object `WaveletCoeffs` is a tuple with ten entries: the first one is the approximation at the highest level: 9; the second entry is a 3-tuple containing the three different details (horizontal, vertical and diagonal) at the same level. Each consecutive entry is a 3-tuple containing the three different details at lower levels

The module `pywt` has several threshold functions implemented. Let us use the soft thresholding indicated above, with a made-up threshold

In order to apply this threshold to each single wavelet coefficient, we will make use of the powerful tools of functional programming from `sage`. We need:

- A thresholding function. We use the ones provided by the
`pywt`package:`pywt.thresholding.soft`, and in particular its lambda expression:lambda x: pywt.thresholding.soft(x, threshold)

`Python`‘s built-in function`map(func, seq1, seq2, ...)`- The tuple of objects over which we need to apply the thresholding: In this case, the set of wavelet coefficients
`WaveletCoeffs`computed above.

Once computed the new coefficients, we use `pywt.waverec2` to obtain the corresponding reconstruction:

threshold = noiseSigma*sqrt(2*log2(image.size)) NewWaveletCoeffs = map (lambda x: pywt.thresholding.soft(x,threshold), WaveletCoeffs) NewImage = pywt.waverec2( NewWaveletCoeffs, wavelet)

original + noise | denoised (haar) |

Denoising is good, but the overall visual quality of the resulting image is very poor: there are many artifacts, due mainly to the blocky nature of the Haar wavelet. Also, the choice of threshold is artificial: it has a priori no relationship with the structure of the image. Let’s put a temporary solution to this problem by using a smoother wavelet, while keeping the same threshold. To avoid excessive typing, we will wrap the code above into a neat function with variables being the image to be treated, the standard deviation of the noise added, and the name of the wavelet employed:

def denoise(data,wavelet,noiseSigma): levels = Integer(floor(log2(data.shape[0]))) WC = pywt.wavedec2(data,wavelet,level=levels) threshold=noiseSigma*sqrt(2*log2(data.size)) NWC = map(lambda x: pywt.thresholding.soft(x,threshold), WC) return pywt.waverec2( NWC, wavelet)

We can now perform the same operation with different wavelets. The following script creates a `python` dictionary that assigns, to each wavelet, the corresponding denoised version of the corrupted Lena image.

Denoised={} for wlt in pywt.wavelist(): Denoised[wlt] = denoise( data=image, wavelet=wlt, noiseSigma=16.0)

The four images below are the respective denoising by soft thresholding of wavelet coefficients on the same image with the same level of noise for the symlet `sym15`, the Daubechies wavelet `db6`, the biorthogonal wavelet `bior2.8`, and the coiflet `coif2`. Note that, except in the case of the denoising by biorthogonal wavelet, all of the others present clear artifacts:

denoised (sym15) |
denoised (db6) |

denoised (bior2.8) |
denoised (coif2) |

This is a great tutorial except that I have no idea what the Integer object is. I’m not sure I need full blown Sage and already have all the other prerequisites installed, including PyWavelets. I’m building sage from source now but that seems a bit crazy to me.

I agree! You definitely do not need the whole

sageto perform numerical methods; just the needed python packages. The nice thing aboutsageis that you can do all sorts ofothermath with it as well. It comes in handy when you need some Graph Theory or big algebraic manipulations (as in manipulations of ideals, for instance). Thanks for your interest!perhaps “levels = (floor(log2(data.shape[0]))).astype(int)” works instead of “levels = Integer(floor(log2(data.shape[0])))”

I still have a question about : threshold=noiseSigma*sqrt(2*log2(data.size))

how to define the noiseSigma provided that the noise level is not clear?

nice but i am working with estimation of wavelet coefficient, so i need a information about that estimation

Great tutorial, very useful and well-explained! Thanks for making it available!!

I would like to implement the mexican hat wavelet in a similar example, but I could not seem to find this particular wavelet available in pywt nor did my queries on the net help me on how to implement it… I thought it was widely used in many branches of science…

Would you have any suggestion on how to do it or where to find a code snippet for that?

Many thanks!!

Thanks for your kind words, Catarina!

It is relatively easy to construct your own wavelet with pywt. All you need to know is their wavelet filter coefficients, for the low- and high-pass decomposition filters. Once you have them, you collect them in an array (let’s call it mexh_filter_bank) and you create the new wavelet by issuing something like

mexh = pywt.wavelet(‘Mexican Hat’, filter_bank=mexh_filter_bank)

(see http://www.pybytes.com/pywavelets/regression/wavelet.html for more detailed info if you need)