Wavelet Denoising¶
This numerical tour uses wavelets to perform non-linear image denoising. $\newcommand{\dotp}[2]{\langle #1, #2 \rangle}$ $\newcommand{\qandq}{\quad\text{and}\quad}$ $\newcommand{\qwhereq}{\quad\text{where}\quad}$ $\newcommand{\qifq}{ \quad \text{if} \quad }$ $\newcommand{\ZZ}{\mathbb{Z}}$ $\newcommand{\RR}{\mathbb{R}}$ $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\si}{\sigma}$ $\newcommand{\Nn}{\mathcal{N}}$ $\newcommand{\Bb}{\mathcal{B}}$ $\newcommand{\EE}{\mathbb{E}}$ $\newcommand{\norm}[1]{\|#1\|}$ $\newcommand{\abs}[1]{\left|#1\right|}$ $\newcommand{\choice}[1]{ \left\{ \begin{array}{l} #1 \end{array} \right. }$ $\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}$
Important: Please read the installation page for details about how to install the toolboxes.
from __future__ import division
from nt_toolbox.general import *
from nt_toolbox.signal import *
%load_ext autoreload
%autoreload 2
%pylab inline
%matplotlib inline
Populating the interactive namespace from numpy and matplotlib
WARNING: pylab import has clobbered these variables: ['pylab'] `%matplotlib` prevents importing * from pylab and numpy
Image Denoising¶
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.
n = 256
name = 'nt_toolbox/data/flowers.png'
f0 = load_image(name, n)
imageplot(f0, 'Original image')
Standard deviation $\si$ of the noise.
sigma = .08
Then we add Gaussian noise $w$ to obtain $f=f_0+w$.
f = f0 + sigma*random.standard_normal(f0.shape)
Display the noisy image.
imageplot(clamp(f), 'Noisy, SNR=' + str(snr(f0,f)) )
Hard Thresholding in Wavelet Bases¶
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$.
def thresh_hard(u,t):return u*(abs(u)>t)
alpha = linspace(-3,3,1000)
plot(alpha, thresh_hard(alpha,1))
axis('equal');
Parameters for the orthogonal wavelet transform.
h = [0, .482962913145, .836516303738, .224143868042, -.129409522551]
h = h/norm(h)
Jmin = 2
First we compute the wavelet coefficients $a$
of the noisy image $f$.
a = perform_wavortho_transf(f,Jmin,+1,h)
Display the noisy coefficients.
plot_wavelet(a,Jmin);
Select the threshold value, that should be proportional to the noise level $\si$.
T = 3*sigma
Hard threshold the coefficients below the noise level to obtain $a_T(m)=s_T^0(a_m)$.
aT = thresh_hard(a,T)
Display the thresholded coefficients.
plot_wavelet(aT,Jmin);
Reconstruct the image $\tilde f$ from these noisy coefficients.
fHard = perform_wavortho_transf(aT,Jmin,-1,h)
Display the denoising result.
imageplot(clamp(fHard), 'Hard, SNR=' + str(snr(f0,fHard)) )
Wavelet Denoising with Soft Thesholding¶
The estimated image $\tilde f$ using hard thresholding. suffers from many artifacts. It is possible to improve the result by using soft thresholding, defined as $$ \tilde f = S_T^1(f) = \sum_m s_T^1(\dotp{f}{\psi_m}) \psi_m $$ $$ \qwhereq s_T^1(\alpha) = \max\pa{0, 1 - \frac{T}{\abs{\alpha}}}\alpha. $$
Display the soft thresholding function $s_T^1(\al)$.
def thresh_soft(u,t):return maximum(1-t/abs(u), 0)*u
alpha = linspace(-3,3,1000)
plot(alpha, thresh_soft(alpha,1))
axis('equal');
Select the threshold.
T = 3/2*sigma
Perform the soft thresholding.
aT = thresh_soft(a,T)
To slightly improve the soft thresholding performance, we do not threshold the coefficients corresponding to coarse scale wavelets.
aT[:2^Jmin:,:2^Jmin:] = a[:2^Jmin:,:2^Jmin:]
Re-construct the soft thresholding estimator $\tilde f$.
fSoft = perform_wavortho_transf(aT,Jmin,-1,h)
Display the soft thresholding denoising result.
imageplot(clamp(fSoft), 'Soft, SNR=' + str(snr(f0,fSoft)) )