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}$
The simplest noise model is Gaussian additive noise, where the variance of the pixel value is independent of the mean (which is the value we look for).
Unfortunately, most real-life data corresponds to noise model that are much more complicated, and in particular the variance of the noise often depends on the parameter of interest (for instance the mean).
from __future__ import division
import nt_toolbox as nt
from nt_solutions import denoisingwav_5_data_dependent as solutions
%matplotlib inline
%load_ext autoreload
%autoreload 2
A Poisson model assume that each pixel $x$ of an image $f(x)$ is drawn from a Poisson distribution of parameter $\lambda=f_0(x)$, where $f_0$ is the clean intensity image to recover.
$$ \PP(f(x)=k)=\lambda^k e^{-\lambda}/k! $$Display the Poisson distribution for several value of $\lambda$.
lambda = [4 10 20]
[k, Lambda] = meshgrid(1: 50, lambda)
P = Lambda.^k .* exp(-Lambda)./ factorial(k)
h = plot(P'); axis('tight')
if using_matlab()
set(h, 'LineWidth', 2)
legend('\lambda = 2', '\lambda = 10', '\lambda = 20')
set_label('k', 'P(k)')
This model corresponds to a photon count, where $\lambda$ is proportional to the number of photons that hits the receptor during the exposition time. This is useful to model medical imaging (confocal microscopy), TEP and SPECT tomography, and digital camera noises.
The goal of denoising is to retrieve the value of $f_0(x)=\lambda(x)$, the true image value, which corresponds to a clean image.
Note that $\lambda$ is the mean of the distribution, but is also its variance, so that the noise intensity perturbating the image pixel $f(x)$ is proportional to $f_0(x)$.
We load a clean, unquantized image.
n = 256
name = 'lena'
f0u = rescale(load_image(name, n))
Quantize the values to given $\lambda$ range.
lmin = 1
lmax = 40
f0 = floor(rescale(f0u, lmin, lmax))
Generate the noisy image.
f = poissrnd(f0)
Display.
imageplot(f0, 'Intensity map f0', 1, 2, 1)
imageplot(f, 'Observed noisy image f', 1, 2, 2)
Display the difference, which shows that the noise level depends on the intensity. The noise is larger in bright areas.
imageplot(f0, 'Intensity map f0', 1, 2, 1)
imageplot(f-f0, 'f-f0', 1, 2, 2)
Exercise 1
Display noisy image contaminated by Poisson noise of varying range.
solutions.exo1()
## Insert your code here.
Parameters for the wavelet transform.
Jmin = 4
options.ti = 1
A translation invariance wavelet transform denoising computes an estimate $\tilde f$ of the clean image $f_0$ as $$ \tilde f = \Psi^+ \circ S_T \circ \Psi f $$ where $\Psi$ is the wavelet transform and $S_T$ the hard thresholding operator using a well chosen threshold $T$.
Exercise 2
Perform translation invariant wavelet hard thresholding directly on the Poisson noisy image $f$. Check for an optimal threshold that maximize the SNR. Record the optimal result in |fPoisson|. isplay est result
solutions.exo2()
## Insert your code here.
Display.
imageplot(f0, 'Original image', 1, 2, 1)
imageplot(clamp(fPoisson, min(f0(: )), max(f0(: ))), strcat(['Denoised, SNR = ' num2str(snr(f0, fPoisson), 4)]), 1, 2, 2)
A variance stabilization transform (VST) is a mapping $\phi$ applied to a noisy image $f$ so that the distribution of each pixel $\phi(f(x))$ is approximately Gaussian.
The Anscombe variance stabilization transform (VST) is given by the non-linear mapping.
$$\phi(a) = 2 \sqrt{a+3/8}$$It transforms a Poisson distribution $P(\lambda)$ to a distribution with approximate variance 1, whatever $\lambda$ is, for $\lambda$ large enough (say $\lambda>20$).
Other VST have been proposed, for instance the Freeman & Tukey VST, given by
$$\phi(a) = \sqrt{a+1}+\sqrt{a}$$Exercise 3
Display the estimated variance of a Poisson distribution for various $\lambda$ (e.g. 10000 realization for each $\lambda$) and display the variance of a stabilized distribution (here the green curve corresponds to 'Freeman & Tukey' and the blue curve to 'Anscombe'. lot bplot(2,1,1); = plot(v); axis('tight'); using_matlab() set(h, 'LineWidth', 2); d tle('Poisson variance'); t_label('\lambda', 'Variance'); bplot(2,1,2);
solutions.exo3()
## Insert your code here.
To perform denosing, one applies the VST, performs wavelet thresholding as if the data was corrupted by an additive Gaussian noise of variance $\sigma=1$, and then applies the inverse VST.
This corresponds to computing an estimate $\tilde f$ of the clean image $f_0$ as $$ \tilde f = \phi^{-1} \pa{ \Psi^+ S_T \circ \Psi \circ \phi(f) } $$ where $\Psi$ is the wavelet transform and $S_T$ the hard thresholding.
Exercise 4
Perform translation invariance wavelet hard thresholding on the variance stabilized image. Use for instance the Anscombe VST. Check for an optimal threshold that maximize the SNR. Record the optimal result in |fVST|. isplay est result
solutions.exo4()
## Insert your code here.
Display.
imageplot(clamp(fPoisson, min(f0(: )), max(f0(: ))), ...
strcat(['Un-stabilized, SNR = ' num2str(snr(f0, fPoisson), 4)]), 1, 2, 1)
imageplot(clamp(fVST, min(f0(: )), max(f0(: ))), ...
strcat(['Stabilized, SNR = ' num2str(snr(f0, fVST), 4)]), 1, 2, 2)
There is no close form solution for the Freeman VST. For each denoised stabilized pixel value $b \in \RR$, one should use a Newton algorithm to solve the equation $\phi(a)=b$ and recovered the denoised value $a \in \RR$. This reads $$ a_{k+1} = a_k - \frac{\phi(a_k)-y}{\phi'(a_k)} $$
Exercise 5
Perform VST denoising using the Freeman VST.
solutions.exo5()
## Insert your code here.
A multiplicative noise corresponds to a noise model where the clean image $f_0$ is multiplied by the noise $W$ to obtained the noisy image $f(x)=W(x) f_0(x)$.
This model is useful for data acquired with an active acquisition device, for instance SAR imaging and ultra-sound imaging.
The distribution of $f(x)$ thus has mean $f_0(x)$ and variance proportional to $f_0(x) \sigma$ where $\sigma$ is the variance of $W$.
Load a clean image.
n = 256
name = 'boat'
f0 = rescale(load_image(name, n), 1e-2, 1)
A classical model for the multiplier noise $W$, that is used for instance in SAR imaging, is to assume that the noise has a Gamma law of mean 1 and variance parameterized by the noise level $L$. $$ P(W=x) \sim x^{K-1} e^{ -x/\theta } $$
where the mean of $P$ is $s=K \theta$ and the variance is $\sigma^2=K \theta^2$.
This corresponds to a acquisition device which is averaging $K$ measures with 1-sided exponential distribution.
K = 4
sigma = 1/ sqrt(K)
Generate the random multiplier.
W = gamrnd(1/ sigma^2, sigma^2, n, n)
Generate the noisy signal.
f = f0.*W
Display.
imageplot(f0, 'Intensity map f0', 1, 2, 1)
imageplot(f, 'Observed noisy image f', 1, 2, 2)
Display the difference, which shows that the noise level depends on the intensity. The noise is larger in bright areas.
imageplot(f0, 'Intensity map f0', 1, 2, 1)
imageplot(f-f0, 'f-f0', 1, 2, 2)
Exercise 6
Generate several noisy images for several noise levels.
solutions.exo6()
## Insert your code here.
Exercise 7
Perform translation invariance wavelet hard thresholding directly on the noisy image $f=f_0 W$. Check for an optimal threshold that maximize the SNR. Record the optimal result in |fMult|. isplay est result
solutions.exo7()
## Insert your code here.
Display.
imageplot(f0, 'Original image', 1, 2, 1)
imageplot(clamp(fMult, min(f0(: )), max(f0(: ))), strcat(['Denoised, SNR = ' num2str(snr(f0, fMult), 4)]), 1, 2, 2)
An approximate variance stabilization transform consist in taking the log. The $\log(f)-a$ image is then equal to $\log(f_0)$ contaminated by $log(W)-a$ which is not too far from a centered Gaussian distribution if the variance of $W$ is not too large.
The value of $a$ should be chosen as the mean value of the random variable $\log(f)$ so that $\log(W)-a$ has zero mean. There exists close form for the value of $a$ as a function of $\sigma=1/\sqrt{K}$, which we use here, where $\psi$ is the polygamma function.
Important: Scilab user should replace |psi| by |dlgamma|.
a = psi(K) - log(K)
Display stabililized image.
imageplot(f, 'f', 1, 2, 1)
imageplot(clamp(log(f)-a, -2, 2), 'log(f)', 1, 2, 2)
Distribution of the noise multiplier and of the log.
subplot(2, 1, 1)
hist(W(: ), 100)
axis('tight')
subplot(2, 1, 2)
hist(log(W(: ))-a, 100)
axis('tight')
Exercise 8
Perform translation invariance wavelet hard thresholding on the variance stabilized image using the log. Check for an optimal threshold that maximize the SNR. Record the optimal result in |fVST|. isplay est result
solutions.exo8()
## Insert your code here.
Display.
imageplot(clamp(fMult, min(f0(: )), max(f0(: ))), strcat(['Un-stabilized, SNR = ' num2str(snr(f0, fMult), 4)]), 1, 2, 1)
imageplot(clamp(fVST, min(f0(: )), max(f0(: ))), strcat(['Stabilized, SNR = ' num2str(snr(f0, fVST), 4)]), 1, 2, 2)