Data Dependent Noise Models
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).
Contents
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).
Installing toolboxes and setting up the path.
You need to download the general purpose toolbox and the signal toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_general/ and toolbox_signal/ in your directory.
For Scilab user: you must replace the Matlab comment '%' by its Scilab counterpart '//'.
Recommandation: You should create a text file named for instance numericaltour.sce (in Scilab) or numericaltour.m (in Matlab) to write all the Scilab/Matlab commands you want to execute. Then, simply run exec('numericaltour.sce'); (in Scilab) or numericaltour; (in Matlab) to run the commands.
Execute this line only if you are using Matlab.
getd = @(p)path(path,p); % scilab users must *not* execute this
Then you can add these toolboxes to the path.
% Add some directories to the path getd('toolbox_signal/'); getd('toolbox_general/');
Poisson Noise
A Poisson model assume that each pixel (i,j) of an image M(i,j) is drawn from a Poisson distribution of parameter lambda=lambda(i,j)
P(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); end 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 M0(i,j)=lambda(i,j), the true image value, which corresponds to a clean image.
Note that lambda is the mean of P, but is also its variance, so that the noise intensity perturbating the image pixel M(i,j) is proportional to M0(i,j).
We load a clean, unquantized image.
n = 256;
name = 'lena';
M0u = rescale( load_image(name,n) );
Quantize the values to given lambda range.
lmin = 1; lmax = 40; M0 = floor( rescale(M0u,lmin,lmax) );
Generate the noisy image.
M = poissrnd(M0);
Display.
clf; imageplot(M0, 'Intensity map M0', 1,2,1); imageplot(M, 'Observed noisy image M', 1,2,2);
Display the difference, which shows that the noise level depends on the intensity. The noise is larger in bright areas.
clf; imageplot(M0, 'Intensity map M0', 1,2,1); imageplot(M-M0, 'M-M0', 1,2,2);
Exercice 1: (the solution is exo1.m) Display noisy image contaminated by Poisson noise of varying range.
exo1;
Parameters for the wavelet transform.
Jmin = 4; options.ti = 1;
Exercice 2: (the solution is exo2.m) Perform translation invariance wavelet hard thresholding directly on the Poisson noisy image M. Check for an optimal threshold that maximize the SNR. Record the optimal result in Mpoisson.
exo2;
Display.
clf; imageplot(M0, 'Original image', 1,2,1); imageplot(clamp(Mpoisson,min(M0(:)),max(M0(:))), strcat(['Denoised, SNR=' num2str(snr(M0,Mpoisson),4)]), 1,2,2);
Variance Stabilization Poisson Denoising
A variance stabilization transform (VST) is a mapping f applied to a noisy image M so that the distribution of each pixel f(M(i,j)) is approximately Gaussian.
The Anscombe variance stabilization transform (VST) is given by the non-linear mapping.
f(x) = 2*sqrt(x+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).
Exercice 3: (the solution is exo3.m) 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.
exo3;
Other VST have been proposed, for instance the Freeman & Tukey VST, given by
f(x) = sqrt(x+1)+sqrt(x)
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.
Exercice 4: (the solution is exo4.m) 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 Mvst.
exo4;
Display.
clf; imageplot(clamp(Mpoisson,min(M0(:)),max(M0(:))), ... strcat(['Un-stabilized, SNR=' num2str(snr(M0,Mpoisson),4)]), 1,2,1); imageplot(clamp(Mvst,min(M0(:)),max(M0(:))), ... strcat(['Stabilized, SNR=' num2str(snr(M0,Mvst),4)]), 1,2,2);
Multiplicative Noise
A multiplicative noise corresponds to a noise model where the clean image M0 is multiplied by the noise W to obtained the noisy image M=W.*M0.
This model is useful for data acquired with an active acquisition device, for instance SAR imaging and ultra-sound imaging.
The distribution of M(i,j) thus has mean M0(i,j) and variance proportional to M0(i,j)*sigma where sigma is the variance of W.
Load a clean image.
n = 256;
name = 'boat';
M0 = 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) ~ x^(k-1)*exp(-x/theta)
where the mean of P is s=k*theta and the variance is sigma^2=k*theta^2.
Variance of the noise multiplier.
sigma = .3;
Generate the random multiplier.
W = gamrnd(1/sigma^2, sigma^2,n,n);
Generate the noisy signal.
M = M0.*W;
Display.
clf; imageplot(M0, 'Intensity map M0', 1,2,1); imageplot(M, 'Observed noisy image M', 1,2,2);
Display the difference, which shows that the noise level depends on the intensity. The noise is larger in bright areas.
clf; imageplot(M0, 'Intensity map M0', 1,2,1); imageplot(M-M0, 'M-M0', 1,2,2);
Exercice 5: (the solution is exo5.m) Generate several noisy images for several noise levels.
exo5;
Exercice 6: (the solution is exo6.m) Perform translation invariance wavelet hard thresholding directly on the noisy image M=M0.*W. Check for an optimal threshold that maximize the SNR. Record the optimal result in Mmult.
exo6;
Display.
clf; imageplot(M0, 'Original image', 1,2,1); imageplot(clamp(Mmult,min(M0(:)),max(M0(:))), strcat(['Denoised, SNR=' num2str(snr(M0,Mmult),4)]), 1,2,2);
Variance Stabilization for Multiplicative Noise
An approximate variance stabilization transform consist in taking the log. The log(M)-a image is then equal to log(M0) 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(M) so that log(W)-a has zero mean. There exists close form for the value of a as a function of sigma, but here we cheat and use directly the estimated mean.
a = mean(log(W(:)));
Display stabililized image.
clf; imageplot(M, 'M', 1,2,1); imageplot(clamp(log(M)-a,-2,2), 'log(M)', 1,2,2);
Distribution of the noise multiplier and of the log.
clf; subplot(2,1,1); hist(W(:),100); axis('tight'); subplot(2,1,2); hist(log(W(:))-a,100); axis('tight');
Exercice 7: (the solution is exo7.m) 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 Mvst.
exo7;
Display.
clf; imageplot(clamp(Mmult,min(M0(:)),max(M0(:))), strcat(['Un-stabilized, SNR=' num2str(snr(M0,Mmult),4)]), 1,2,1); imageplot(clamp(Mvst,min(M0(:)),max(M0(:))), strcat(['Stabilized, SNR=' num2str(snr(M0,Mvst),4)]), 1,2,2);
