This numerical tour introduces basic image denoising methods.
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}$
from __future__ import division
from nt_toolbox.general import *
from nt_toolbox.signal import *
%pylab inline
%matplotlib inline
%load_ext autoreload
%autoreload 2
/Users/gpeyre/anaconda/lib/python3.5/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment. warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.') /Users/gpeyre/anaconda/lib/python3.5/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment. warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')
Populating the interactive namespace from numpy and matplotlib
/Users/gpeyre/anaconda/lib/python3.5/site-packages/IPython/core/magics/pylab.py:161: UserWarning: pylab import has clobbered these variables: ['pylab'] `%matplotlib` prevents importing * from pylab and numpy "\n`%matplotlib` prevents importing * from pylab and numpy"
In these numerical tour, we simulate noisy acquisition by adding some white noise (each pixel is corrupted by adding an independant Gaussian variable).
This is useful to test in an oracle maner the performance of our methods.
Size $N = n \times n$ of the image.
n = 256
N = n**2
We load a clean image $x_0 \in \RR^N$.
name = 'nt_toolbox/data/flowers.png'
x0 = load_image(name, n)
Display the clean image.
imageplot(x0)
Variance of the noise.
sigma = .08
We add some noise to it to obtain the noisy signal $y = x_0 + w$. Here $w$ is a realization of a Gaussian white noise of variance $\si^2$.
y = x0 + sigma*random.standard_normal(x0.shape)
Display the noisy image.
imageplot(clamp(y))
We consider a noising estimator $x \in \RR^N$ of $x_0$ that only depends on the observation $y$. Mathematically speaking, it is thus a random vector that depends on the noise $w$.
A translation invariant linear denoising is necessarely a convolution with a kernel $h$ $$ x = x_0 \star h $$ where the periodic convolution between two 2-D arrays is defined as $$ (a \star b)_i = \sum_j a(j) b(i-j). $$
It can be computed over the Fourier domain as $$ \forall \om, \quad \hat x(\om) = \hat x_0(\om) \hat h(\om). $$
cconv = lambda a,b : real(ifft2(fft2(a)*fft2(b)))
We use here a Gaussian fitler $h$ parameterized by the bandwith $\mu$.
normalize = lambda h : h/sum(h.flatten())
t = transpose( concatenate( (arange(0,n/2), arange(-n/2,0) ) ) )
[Y,X] = meshgrid(t, t)
h = lambda mu: normalize(exp(-(X**2 + Y**2)/ (2*mu**2)))
Display the filter $h$ and its Fourier transform.
mu = 10
subplot(1,2, 1)
imageplot(fftshift(h(mu)))
title('h')
subplot(1,2, 2)
imageplot(fftshift(real(fft2(h(mu)))))
title('$\hat h$');
imageplot(h(mu))
Shortcut for the convolution with $h$.
denoise = lambda x,mu : cconv(h(mu), x)
Display a denoised signal.
imageplot(denoise(y, mu))
Exercise 1: Display a denoised signal for several values of $\mu$.
run -i nt_solutions/denoisingsimp_2b_linear_image/exo1