Advanced Wavelet Thresholdings

Important: Please read the installation page for details about how to install the toolboxes. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\enscond}[2]{\lbrace #1, #2 \rbrace}$ $\newcommand{\pd}[2]{ \frac{ \partial #1}{\partial #2} }$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\umax}[1]{\underset{#1}{\max}\;}$ $\newcommand{\umin}[1]{\underset{#1}{\min}\;}$ $\newcommand{\uargmin}[1]{\underset{#1}{argmin}\;}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\diag}[1]{{diag}\left( #1 \right)}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\qarrq}{ \quad \Longrightarrow \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\CC}{\mathbb{C}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\Zz}{\mathcal{Z}}$ $\newcommand{\Ww}{\mathcal{W}}$ $\newcommand{\Vv}{\mathcal{V}}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\NN}{\mathcal{N}}$ $\newcommand{\Hh}{\mathcal{H}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\Ee}{\mathcal{E}}$ $\newcommand{\Cc}{\mathcal{C}}$ $\newcommand{\Gg}{\mathcal{G}}$ $\newcommand{\Ss}{\mathcal{S}}$ $\newcommand{\Pp}{\mathcal{P}}$ $\newcommand{\Ff}{\mathcal{F}}$ $\newcommand{\Xx}{\mathcal{X}}$ $\newcommand{\Mm}{\mathcal{M}}$ $\newcommand{\Ii}{\mathcal{I}}$ $\newcommand{\Dd}{\mathcal{D}}$ $\newcommand{\Ll}{\mathcal{L}}$ $\newcommand{\Tt}{\mathcal{T}}$ $\newcommand{\si}{\sigma}$ $\newcommand{\al}{\alpha}$ $\newcommand{\la}{\lambda}$ $\newcommand{\ga}{\gamma}$ $\newcommand{\Ga}{\Gamma}$ $\newcommand{\La}{\Lambda}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Si}{\Sigma}$ $\newcommand{\be}{\beta}$ $\newcommand{\de}{\delta}$ $\newcommand{\De}{\Delta}$ $\newcommand{\phi}{\varphi}$ $\newcommand{\th}{\theta}$ $\newcommand{\om}{\omega}$ $\newcommand{\Om}{\Omega}$

This numerical tour present some advanced method for denoising that makes use of some exoting 1D thresholding functions, that in some cases give better results than soft or hard thresholding.

In [2]:
warning off
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/denoisingwav_3_advanced')
warning on

Generating a Noisy Image

Here we use an additive Gaussian noise.

First we load an image.

In [3]:
n = 256;
name = 'hibiscus';
M0 = load_image(name,n);
M0 = rescale( sum(M0,3) );

Noise level.

In [4]:
sigma = .08;

Then we add some Gaussian noise to it.

In [5]:
M = M0 + sigma*randn(size(M0));

Compute a 2D orthogonal wavelet transform.

In [6]:
Jmin = 3;
MW = perform_wavelet_transf(M,Jmin,+1);

Semi-soft Thresholding

Hard and soft thresholding are two specific non-linear diagonal estimator, but one can optimize the non-linearity to capture the distribution of wavelet coefficient of a class of images.

Semi-soft thresholding is a familly of non-linearities that interpolates between soft and hard thresholding. It uses both a main threshold |T| and a secondary threshold |T1=mu*T|. When |mu=1|, the semi-soft thresholding performs a hard thresholding, whereas when |mu=infty|, it performs a soft thresholding.

In [7]:
T = 1; % threshold value
v = linspace(-5,5,1024);
clf;
hold('on');
plot(v, perform_thresholding(v,T,'hard'), 'b--');
plot(v, perform_thresholding(v,T,'soft'), 'r--');
plot(v, perform_thresholding(v,[T 2*T],'semisoft'), 'g');
plot(v, perform_thresholding(v,[T 4*T],'semisoft'), 'g:');
legend('hard', 'soft', 'semisoft, \mu=2', 'semisoft, \mu=4');
hold('off');

Exercise 1

Compute the denoising SNR for different values of |mu| and different value of |T|. Important: to get good results, you should not threshold the low frequency residual. list/sigma, mulist,

In [8]:
exo1()
In [9]:
%% Insert your code here.

One can display, for each |mu|, the optimal SNR (obtained by testing many different |T|). For |mu=1|, one has the hard thresholding, which gives the worse SNR. The optimal SNR is atained here for |mu| approximately equal to 6.

In [10]:
err_mu = compute_min(err, 2);
clf;
plot(mulist, err_mu, '.-');
axis('tight');
set_label('\mu', 'SNR');

Soft and Stein Thresholding

Another way to achieve a tradeoff between hard and soft thresholding is to use a soft-squared thresholding non-linearity, also named a Stein estimator.

We compute the thresholding curves.

In [11]:
T = 1; % threshold value
v = linspace(-4,4,1024);

hard thresholding

In [12]:
v_hard = v .* (abs(v)>T);

soft thresholding

In [13]:
v_soft = v .* max( 1-T./abs(v), 0 );

Stein thresholding

In [14]:
v_stein = v .* max( 1-(T^2)./(v.^2), 0 );

We display the classical soft/hard thresholders.

In [15]:
clf;
hold('on');
plot(v, v_hard, 'b');
plot(v, v_soft, 'r');
plot(v, v_stein, 'k--');
hold('off');
legend('Hard', 'Soft', 'Stein');

Exercise 2

Compare the performance of Soft and Stein thresholders, by determining the best threshold value.

In [16]:
exo2()
In [17]:
%% Insert your code here.