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 overviews the use of wavelets for image approximation and denoising.
from __future__ import division
import nt_toolbox as nt
from nt_solutions import introduction_5_wavelets_2d as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
Image approximation is obtained by thresholding wavelets coefficients.
First we load an image $f \in \mathbb{R}^N$ of $N = n_0 \times n_0$ pixels.
name = 'cortex'
n0 = 512
f = load_image(name, n0)
f = rescale(sum(f, 3))
Display it.
imageplot(f)
An orthogonal wavelet basis $ \mathcal{B} = \{ \psi_{j,n}^k \}_{j,n} $ of $\mathbb{R}^N$ is composed of multiscale atoms parameterized by their scale $2^j$, position $2^j n \in [0,1]^2$ and orentation $ k \in \{H,V,D\}$.
A forward wavelet transform computes the set of inner products $$ \Psi f = \{ \langle f,\psi_{j,n}^k\rangle \} \in \mathbb{R}^N $$ where the wavelet atoms are defined as $$ \psi_{j,n}^k(x) = \psi^k\left( \frac{x - 2^j n}{2^j} \right). $$
Set the minimum scale for the transform.
Jmin = 0
A short-cut for the wavelet transform $\Psi$:
Psi = lambda f: perform_wavelet_transf(f, Jmin, + 1)
A short-cut for the inverse wavelet transform $\Psi^{-1} = \Psi^*$:
PsiS = lambda fw: perform_wavelet_transf(fw, Jmin, -1)
Perform the wavelet transform to compute $\Psi f$.
fW = Psi(f)
Display the transformed coefficients.
plot_wavelet(fW)
To perform non-linear image approximation, one remove the small amplitude coefficients. This is performed using a hard thresholding $$ H_T(f,\mathcal{B}) = \Psi^{-1} \circ H_T \circ \Psi (f) = \sum_{|\langle f,\psi_{j,n}^k\rangle| > T} \langle f,\psi_{j,n}^k\rangle \psi_{j,n}^k. $$
$T$ should be adapted to ensure a given number $M$ of non-zero coefficients, and then $f_M = H_T(f,\mathcal{B})$ is the best $M$ terms approximation of $f$ in $\mathcal{B}$.
Select a threshold.
T = .5
Shortcut for the thresholding operator $H_T$.
Thresh = lambda fW, T: fW .* (abs(fW) >T)
Perform hard thresholding of the coefficients.
fWT = Thresh(fW, T)
Exercise 1
Compute the ratio $M/N$ of non-zero coefficients.
solutions.exo1()
## Insert your code here.
Display the thresholded coefficients.
subplot(1, 2, 1)
plot_wavelet(fW)
title('Original coefficients')
subplot(1, 2, 2)
plot_wavelet(fWT)
Perform reconstruction using the inverse wavelet transform $\Psi^*$.
f1 = PsiS(fWT)
Display approximation.
imageplot(f, 'Image', 1, 2, 1)
imageplot(clamp(f1), strcat(['Approximation, SNR = ' num2str(snr(f, f1), 3) 'dB']), 1, 2, 2)
Number of coefficients for the approximation.
M = n0^2/ 16
Exercise 2
Compute a threshold $T$ to keep only $M$ coefficients.
solutions.exo2()
## Insert your code here.
Perform hard thresholding.
fWT = Thresh(fW, T)
Check the number of non-zero coefficients in |fWT|.
disp(strcat([' M = ' num2str(M)]))
disp(strcat(['|fWT|_0 = ' num2str(sum(fWT(: )~ = 0))]))
Exercise 3
Compute an approximation with an decreasing number of coefficients.
solutions.exo3()
## Insert your code here.
Image denoising is obtained by thresholding noisy wavelets coefficielts.
Here we consider a simple setting where we intentionnaly add some noise $w$ to a clean image $f$ to obtain $ y = f + w $.
A denoiser computes an estimate $\tilde f$ of $f$ from the observations $y$ alone. In the mathematical model, since $y$ is a random variable depending on $w$, so is $\tilde f$. A mathematical evaluation of the efficiency of the denoiser is the average risk $E_w( \|f-\tilde f\|^2 )$.
Here we consider a single realization of the noise, so we replace the risk by the oracle error $ \|f-\tilde f\|^2$. This allows us to bench the efficiency of the denoising methods by comparing the result to $f$. But you have to keep in mind that for real application, one does not have access to $f$.
We consider a Gaussian white noise $w$ of variance $\sigma^2$.
sigma = .1
Generate a noisy image.
y = f + randn(n0, n0)*sigma
Display.
imageplot(f, 'Clean image', 1, 2, 1)
imageplot(clamp(y), ['Noisy image, SNR = ' num2str(snr(f, y), 3) 'dB'], 1, 2, 2)
A denoising is obtained by thresholding the wavelet coefficients $$ \tilde f = H_T(f,\mathcal{B}). $$
The asymptotically optimal threshold of Donoho and Johnstone is $T = \sqrt{2 \log(N)} \sigma$. In practice, one observes that much better result are obtained using $T \approx 3 \sigma$.
Compute the noisy wavelet coefficients.
fW = Psi(y)
Compute the threshold value using the $3\sigma$ heuristic.
T = 3*sigma
Perform hard thresholding.
fWT = Thresh(fW, T)
Display the thresholded coefficients.
subplot(1, 2, 1)
plot_wavelet(fW)
title('Original coefficients')
subplot(1, 2, 2)
plot_wavelet(fWT)
Perform reconstruction.
f1 = PsiS(fWT)
Display denoising.
imageplot(clamp(y), 'Noisy image', 1, 2, 1)
imageplot(clamp(f1), strcat(['Denoising, SNR = ' num2str(snr(f, f1), 3) 'dB']), 1, 2, 2)
Exercise 4
Try to optimize the value of the threshold $T$ to get the best possible denoising result.
solutions.exo4()
## Insert your code here.
The quality of orthogonal denoising is improved by adding translation invariance. This corresponds to denoising translated copies of the image.
The translation of an image is $(\theta_\tau f)(x) = f(x-\tau)$, where we use periodic boundary conditions.
Given a set $ \Omega \subset \mathbb{R}^2 $, the $\Omega$-translation invariant denoising is defined as: $$ \tilde f = \frac{1}{\Omega}\sum_{\tau \in \Omega} \theta_{-\tau} \left( H_T( \theta_\tau y, \mathcal{B} ) \right). $$
Here we consider translation of integer pixels in $\{0,\ldots,\tau_{\max}-1\}$. The number of translations is thus $ \tau_{\max}^2$.
tau_max = 8
Generate a set of translation vectors $\Omega = \{ \tau_i = (X_i,Y_i) \}_i$.
[Y, X] = meshgrid(0: tau_max-1, 0: tau_max-1)
A "trick" to compute the full denoising image after all translations is to initialize $\tilde f = 0$, and then accumulate each denoising with translate $\tau_i$ in the following way: $$ \tilde f \longleftarrow \frac{i-1}{i} \tilde f + \frac{1}{i} \theta_{-\tau_i} \left( H_T( \theta_{\tau_i} y, \mathcal{B} ) \right) $$
Initialize the denoised image $\tilde f$ as 0.
f1 = zeros(n0, n0)
Initialize the translation index.
i = 1
Translate the image to obtain $\theta_{\tau_i}(f)$ for $\tau_i = (X_i,Y_i)$, with periodic boundary conditions.
fTrans = circshift(y, [X(i) Y(i)])
Denoise this translated image, to obtain $H_T(\theta_{\tau_i} f,\mathcal{B})$.
fTrans = PsiS(Thresh(Psi(fTrans) , T))
Translate back.
fTrans = circshift(fTrans, -[X(i) Y(i)])
Accumulate the result.
f1 = (i-1)/ i*f1 + fTrans/ i
Exercise 5
Compute the full denoising by cycling through the $i$ indices.
solutions.exo5()
## Insert your code here.
Exercise 6
Determine the optimal threshold $T$ for this translation invariant denoising.
solutions.exo6()
## Insert your code here.
Exercise 7
Test on other images.
solutions.exo7()
## Insert your code here.