Outliers and Median Denoiser
This numerical tour explores non-linear denoisers based on local median computations.
Contents
Installing toolboxes and setting up the path.
You need to download the following files: signal toolbox and general toolbox.
You need to unzip these toolboxes in your working directory, so that you have toolbox_signal and toolbox_general 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 command 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(p,path); % scilab users must *not* execute this
Then you can add the toolboxes to the path.
getd('toolbox_signal/'); getd('toolbox_general/');
Robust Statistics and Median Estimator
Robust statistic are estimators that can compute accurately means and standard deviations of process contaminated with outliers (few exceptional events like spikes). This is related to replacing the L2 norm for which traditional estimator (empirical mean and empirical deviation) are optimal by L1 norm which is robust to outliers.
First we generate a white noise of mean and deviation 1 with a few outliers that have a much larger deviation and that are positive.
% dimension n = 10000; % deviation for base signal and outliers sigma0 = 1; sigma1 = 20; % ratio of outliers rho = .02; % number of outliers p = round(rho*n); % noisy base signal x = randn(n,1)*sigma0; % positions of the outliers sel = randperm(n); sel = sel(1:p); % add outliers x(sel) = x(sel) + rand(p,1)*sigma1; % display clf; plot(x); axis('tight');
We can compute the mean and the median. Note how the mean is shifted because of the positive outliers.
disp(strcat(['Mean =' num2str(mean(x))] )); disp(strcat(['Median=' num2str(median(x))] ));
Mean =0.20045 Median=0.03191
We can compute the standard deviation and the L1 norm of absolute deviation (mad). The mad estimator is normalized by
disp(strcat(['Mean =' num2str(std(x))] )); disp(strcat(['MAD =' num2str(mad(x,1)/0.6745)] ));
Mean =1.8631 MAD =1.0254
Another options to stabilize the estimation is to threshold the signal to detect outliers. In wavelet denoising, such a technic is used because the outlier are in fact the data one is interested in (the signal content). The Gaussian variable is with high probabilitiy below T=sqrt(2*log(n)), so that this universal threshold can be used. In practice, a smaller, less conservative threshold can be used.
T = sqrt(2*log(n)); xT = x .* (abs(x)<=T); plot(xT); axis('tight'); title('Signal below the threshold'); % compute the thresholded mean and deviation disp(strcat(['Mean =' num2str(mean(xT))] )); disp(strcat(['Std =' num2str(std(xT))] ));
Mean =0.016521 Std =0.99715
Signal Impulse Noise Removal with Local Median
A median is able to corretly estimate the mean of a signal. The mean of a smooth signal is a good approximation of the signal itself. For a piecewise smooth signal, this mean needs to be computed locally. Applying local mean (average) leads to blurring and also does not works to remove impulse noise. Both issues are solved by computing a local median.
We load a piecewise smooth signal.
n = 1024;
name = 'piece-regular';
f0 = load_signal(name, n);
f0 = rescale(f0);
We add an impulse noise to f0.
% ratio of perturbated values rho = .3; % amplitude of the perturbation sigma = 1; % number of perturbated values p = round(rho*n); % indexes of the perturbated values sel = randperm(n); sel = sel(1:p); sel = sel(:); % pertubation of the signal f = f0; f(sel) = f(sel) + randn(size(sel)) * sigma;
We display both signals.
clf; subplot(2,1,1); plot(f0); axis('tight'); title('Original signal'); subplot(2,1,2); plot(clamp(f)); axis('tight'); title('Noisy signal');
The only parameter of the algorithm is the half-width k of the window.
% half width k = 3; % width w = 2*k+1;
We initialize the denoising algorithm using the noisy signal.
fmed = x;
For a given position x in the signal, we compute a window of indexes of size 2*w+1 around x.
% central position x = 100; % posisition around x sel = x-k:x+k; % cyclic boundary conditions sel = mod(sel-1,n)+1;
The denoised signal at location x is obtained as the median.
fmed(x) = median(f(sel));
Exercice 1: (the solution is exo1.m) Implement the denoising algorithm by looping through the position x in the signal.
exo1;
Exercice 2: (the solution is exo2.m) Compute the SNR for various window width w, find the optimal window with. Try for several noise level rho and sigma to see how this optimal width change.
exo2;
Image Impulse Noise Removal with Local Median
The same local median estimator is applied to denoise a 2D image.
Load a clean image.
n = 128;
name = 'hibiscus';
M0 = load_image(name,n*2);
M0 = rescale( sum(M0,3) );
M0 = rescale(crop(M0,n));
Parameters for the noise
rho = .4; sigma = 2;
Exercice 3: (the solution is exo3.m) Compute an image M corrupted by adding impulse gaussian noise of variance sigma to M0 for a fraction rho of its pixels. For the display you need to clamp the images.
exo3;
You can try a Gaussian blur to denoise the image. This is not very efficient because a blur is similar to a local mean that does not corretely estimate the true signal.
s = 5; Mh = perform_blurring(M, s); clf; imageplot(clamp(M), 'Noisy', 1,2,1); imageplot(clamp(Mh), 'Blurred', 1,2,2);
For Scilab, one might need to extend a little the size of the stack. Warning: do this only once.
extend_stack_size(4);
A local median applied over a window of w x w pixels removes the outlier of the image. The parameter w should be large enough for the median to be well estimated (the more noise, the larger w) and small enough not too blur too much the edges.
k = 2; % half window w = 2*k+1; % window size Mmed = perform_median_filtering(M,k); clf; imageplot(clamp(M), 'Noisy', 1,2,1); imageplot(clamp(Mmed), 'Median filtered', 1,2,2);
Exercice 4: (the solution is exo4.m) Try for various patch size w to find the best result. Compare the optimal result Mmed with Mh.
exo4;
Iterated median.
Instead of applying the median only once, one can iteratively filter an image to simulate a time evolution. One can prove that asymptotically, this flow converges (up to renormalizaiton) to the Total Variation (TV) flow of Rudin/Osher/Fatemi.
Exercice 5: (the solution is exo5.m) Implemented an iterated median filtering, starting from the image M0. What do you notice. Compare this median with a total variation flow. Use a small w=3.
exo5;
Exercice 6: (the solution is exo6.m) Display the evolution of the denoising error SNR(M0,Mmed) through the iterations.
exo6;
Color Image Denoising
The numerical tour Color Image Denoising with Median Filtering extend the definition of the median to multidimensional points, and use it to perform color image denoising.