*Important:* Please read the installation page for details about how to install the toolboxes.

In [2]:

```
using PyPlot
using NtToolBox
# using Autoreload
# arequire("NtToolBox")
```

We consider a simple generative model of noisy images $F = f_0+W$ where $f_0 \in \RR^N$ is a deterministic image of $N$ pixels, and $W$ is a Gaussian white noise distributed according to $\Nn(0,\si^2 \text{Id}_N)$, where $\si^2$ is the variance of noise.

The goal of denoising is to define an estimator $\tilde F$ of $f_0$ that depends only on $F$, i.e. $\tilde F = \phi(F)$ where $\phi : \RR^N \rightarrow \RR^N$ is a potentially non-linear mapping.

Note that while $f_0$ is a deterministic image, both $F$ and $\tilde F$ are random variables (hence the capital letters).

The goal of denoising is to reduce as much as possible the denoising error given some prior knowledge on the (unknown) image $f_0$. A mathematical way to measure this error is to bound the quadratic risk $\EE_w(\norm{\tilde F - f_0}^2)$, where the expectation is computed with respect to the distribution of the noise $W$.

Image loading and adding Gaussian Noise For real life applications, one does not have access to the underlying image $f_0$. In this tour, we however assume that $f_0$ is known, and $f = f_0 + w\in \RR^N$ is generated using a single realization of the noise $w$ that is drawn from $W$. We define the estimated deterministic image as $\tilde f = \phi(f)$ which is a realization of the random vector $\tilde F$.

First we load an image $f_0 \in \RR^N$ where $N=n \times n$ is the number of pixels, and then display it.

In [2]:

```
n = 256
name = "NtToolBox/src/data/flowers.png"
f0 = load_image(name, n)
imageplot(f0, "Original image")
```

libpng warning: iCCP: known incorrect sRGB profile libpng warning: iCCP: known incorrect sRGB profile libpng warning: iCCP: known incorrect sRGB profile

Out[2]:

PyObject <matplotlib.text.Text object at 0x329420a10>

Standard deviation $\si$ of the noise.

In [3]:

```
sigma = .08;
```

Then we add Gaussian noise $w$ to obtain $f=f_0+w$.

In [4]:

```
using Distributions
f = f0 + sigma.*rand(Normal(), size(f0)[1], size(f0)[2]);
```

Display the noisy image.

In [5]:

```
imageplot(clamP(f), string("Noisy, SNR = ", string(snr(f0,f)) ))
```

Out[5]:

PyObject <matplotlib.text.Text object at 0x329f71210>

A simple but efficient non-linear denoising estimator is obtained by thresholding the coefficients of $f$ in a well chosen orthogonal basis $\Bb = \{\psi_m\}_m$ of $\RR^N$.

In the following, we will focuss on a wavelet basis, which is efficient to denoise piecewise regular images.

The hard thresholding operator with threshold $T \geq 0$ applied to some image $f$ is defined as $$ S_T^0(f) = \sum_{\abs{\dotp{f}{\psi_m}}>T} \dotp{f}{\psi_m} \psi_m = \sum_m s_T^0(\dotp{f}{\psi_m}) \psi_m $$ where the hard thresholding operator is $$ s_T^0(\alpha) = \choice{ \alpha \qifq \abs{\al}>T, \\ 0 \quad\text{otherwise}. }$ $$

The denoising estimator is then defined as $$ \tilde f = S_T^0(f). $$

Display the function $s_T^0(\al)$ for $T=1$.

In [6]:

```
thresh_hard = (u, t) -> u.*(abs(u) .> t)
alpha = linspace(-3, 3, 1000)
plot(alpha, thresh_hard(alpha, 1))
axis("equal")
```

Out[6]:

(-3.3,3.3,-3.3,3.3)

Parameters for the orthogonal wavelet transform.

In [7]:

```
h = [0, .482962913145, .836516303738, .224143868042, -.129409522551]
h = h./norm(h)
Jmin = 2;
```

First we compute the wavelet coefficients $a$

of the noisy image $f$.

In [8]:

```
a = perform_wavortho_transf(f, Jmin, +1, h);
```

Display the noisy coefficients.

In [9]:

```
plot_wavelet(a, Jmin);
```