Curvelet Denoising¶

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 explores the use of curvelets to perform image denoising.

In [22]:

warning off
addpath('toolbox_signal')
addpath('toolbox_general')
addpath('solutions/denoisingwav_6_curvelets')
warning on

Curvelet Transform¶

The curvelet tight frame was introduced by Candes and Donoho to enhance the wavelet representation of geometric cartoon images. See for instance

E.J. Cand s and D. L. Donoho, New tight frames of curvelets and optimal representations of objects with piecewise-C2 singularities, Comm. Pure Appl. Math., 57 219-266

In this tour, we use the discrete curvelet transform, detailed in

E. J. Cand s, L. Demanet, D. L. Donoho and L. Ying, Fast discrete curvelet transforms, Multiscale Model. Simul., 5 861-899

and it uses the <http://www.curvelet.org/ Curvelab> implementation.

Load an image.

In [3]:

n = 256;
name = 'lena';
M = rescale(load_image(name, n));

Parameters for the curvelet transform.

In [4]:

options.null = 0;
options.finest = 1;
options.nbscales = 4;
options.nbangles_coarse = 16;
options.is_real = 1;
options.n = n;

Perform the transform.

In [5]:

MW = perform_curvelet_transform(M, options);

Display the transform.

In [6]:

clf;
plot_curvelet(MW, options);

One can perform a non-linear approximation of the image by thresholding the curvelet coefficients.

In [7]:

T = .2;
MWT = perform_thresholding(MW, T, 'hard');

One recovers the approximated image by using the inverse curvelet transform.

In [8]:

M1 = perform_curvelet_transform(MWT, options);

Display the approximation. Note that is not the best non-linear approximation since the curvelet tight frame is not an orthogonal basis.

In [9]:

clf;
imageplot(M, 'Original', 1,2,1);
imageplot(clamp(M1), 'Approximated', 1,2,2);

Denoising with the Curvelet transform¶

Curvelet thresholding is useful to perform denoising.

We load a sub-set of the lena image.

In [10]:

name = 'lena';
n = 128;
M0 = rescale(crop(load_image(name),n, [108 200]));
options.n = n;

Add some noise.

In [11]:

sigma = .05;
M = M0 + sigma*randn(n);

Exercise 1

Compute the best threshold to minimize the denoising error in curvelets. Call |Mcurv| the optimal denoising.

In [12]:

exo1()

In [13]:

%% Insert your code here.

Display.

In [14]:

clf;
imageplot(clamp(M), 'Noisy', 1,2,1);
imageplot(clamp(Mcurv), ['Denoised, SNR=' num2str(snr(M0,Mcurv),3) 'dB'], 1,2,2);

Exercise 2

Perform cycle spinning to enhance the recovery error.

In [15]:

exo2()

In [16]:

%% Insert your code here.

Display.

In [17]:

clf;
imageplot(clamp(M), 'Noisy', 1,2,1);
imageplot(clamp(Mcurv), ['Denoised, SNR=' num2str(snr(M0,Mcurv),3)], 1,2,2);

Exercise 3

Compare with translation invariant hard thresholding.

In [18]:

exo3()