Non Local Means
This numerical tour study image denoising using non-local means.
Contents
The NL-means algorithm has been introduced for denoising purposes in
A. Buades, B. Coll, J.M Morel, A review of image denoising algorithms, with a new one, SIAM Multiscale Modeling and Simulation, Vol 4 (2), pp: 490-530, 2005.
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 Scilabe) or numericaltour.m 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(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/');
Patches in Images
This numerical tour is dedicated to the study of the structure of patches in images.
We load a noisy image.
n = 128; M0 = load_image('lena'); c = [100 200]; M0 = rescale( crop(M0,n, c) ); % noise level sigma = .04; % add some noise M = M0 + randn(n,n)*sigma;
Display clean and noisy images (use clamp to avoid contrast rescaling during the display).
clf; imageplot(M0, 'Original image', 1,2,1); imageplot(clamp(M), 'Noisy image', 1,2,2);
We set w the (half) size of the patches.
w = 4; w1 = 2*w+1;
We set up larges (n,n,w1,w1) matrices representing the X and Y position of the pixel to extract.
% location of pixels [Y,X] = meshgrid(1:n,1:n); % offsets [dY,dX] = meshgrid(-w:w,-w:w); % location of pixels to extract dX = reshape(dX, [1 1 w1 w1]); dY = reshape(dY, [1 1 w1 w1]); X = repmat(X, [1 1 w1 w1]) + repmat(dX, [n n 1 1]); Y = repmat(Y, [1 1 w1 w1]) + repmat(dY, [n n 1 1]);
We handle boundary condition by reflexion
X(X<1) = 2-X(X<1); Y(Y<1) = 2-Y(Y<1); X(X>n) = 2*n-X(X>n); Y(Y>n) = 2*n-Y(Y>n);
We extract all the patche in the image.
P = M(X + (Y-1)*n);
Each P(i,j,:,:,) represent an (w1,w1) patch extracted around pixel (i,j) in the image.
Display some example of patches
% the image clf; imageplot(M); % some randomly selected patches clf; for i=1:16 x = floor( rand*(n-1)+1 ); y = floor( rand*(n-1)+1 ); imageplot( squeeze(P(x,y,:,:)), '', 4,4,i ); end
Dimensionality Reduction with PCA
Since NL-means type algorithm require the computation of many distances between patches, it is advantagous to reduce the dimensionality of the patch while keeping as much as possible of information.
Target dimensionality
d = 16;
A linear dimensionality reduction is obtained by Principal Component Analysis (PCA) that projects the data on a small number of leading direction of the coveriance matrix of the patches.
Turn the patch matrix into an (w1*w1,n*n) array, so that each P(:,i) is a w1*w1 vector representing a patch.
P = reshape(P, [n*n w1*w1])';
Compute the mean and the covariance of the points cloud representing the patches.
mu = mean(P); P1 = P - repmat(mu, [w1*w1 1]); C = P1*P1';
Extract the eigenvectors.
[V,D] = eig(C); D = diag(D); % sort by decreasing amplitude [D,I] = sort(D, 'descend'); V = V(:,I);
Display the decaying amplitude of the eigenvalues.
clf;
plot(D); axis('tight');
Display the leading eigenvectors - they look like Fourier modes.
clf; for i=1:16 imageplot( reshape(V(:,i),[w1 w1]), '', 4,4,i ); end
Perform the dimension reduction.
H = V(:,1:d)' * P1;
Reshape the array so that each H(i,j,:) is a d dimensional descriptor of a patch.
H = reshape(H', [n n d]);
Non-local Filter
NL-means applies, to each pixel location, an adaptive averaging kernel that is computed from patch distances.
Position of a pixel.
x = 83; y = 72;
Square distance between the patch around (x,y) and all the other ones.
D = sum( (H - repmat(H(x,y,:), [n n 1])).^2, 3 )/(w1*w1);
The non-local filter is obtained from D using a Gaussian function. The width of this Gaussian is very important and should be adapted to match the noise level.
tau = .05;
Compute and normalize the weight.
K = exp( -D/(2*tau^2) ); K = K/sum(K(:));
Display the distance and the kernel.
clf; imageplot(sqrt(D), 'Distance', 1,2,1); imageplot(K, 'Kernel', 1,2,2);
The NL-filtered value at pixel (x,y) is obtained by averaging the values of M with the weight K.
NLval = sum(sum(K.*M));
Localizing the Non-local Means
Non-local Means computes an average value for all pixels (x,y).
We set a "locality constant" that set the maximum distance between patches to compare. This allows to speed up computation, and makes NL-means type methods semi-global (to avoid searching in all the image).
q = 14;
Using this locality constant, we compute the distance between patches only within a window.
selx = x-q:x+q; sely = y-q:y+q;
Once again, one should be careful about boundary conditions.
selx(selx<1 | selx>n) = []; sely(sely<1 | sely>n) = [];
Compute distance and kernel only within the window.
D = sum( (H(selx,sely,:) - repmat(H(x,y,:), [length(selx) length(sely) 1])).^2, 3 )/(w1*w1); K = exp( -D/(2*tau^2) ); K = K/sum(K(:));
Display the distance and the kernel.
clf; imageplot(M(selx,sely), 'Image', 1,2,1); imageplot(K, 'Kernel', 1,2,2);
The NL-filtered value at pixel (x,y) is obtained by averaging the values of M with the weight K.
NLval = sum(sum(K.*M(selx,sely)));
Non-local Means for Denoising
Non-local means is typically applied to denoise an image.
Exercice 1: (the solution is exo1.m) Compute the Non-local Means denoising Mnl by applying the localized NL filter to each pixel. You can use for instance tau=.02
exo1;
Exercice 2: (the solution is exo2.m) Compute the denoising result for several values of tau in order to determine the optimal denoising Mnl0 that maximizes snr(M0,Mnl).
exo2;
Display optimally denoised result.
clf; imageplot(clamp(M), ['Noisy, SNR=' num2str(snr(M0,M)) 'dB'], 1,2,1); imageplot(clamp(Mnl0), ['Denoised, SNR=' num2str(snr(M0,Mnl0)) 'dB'], 1,2,2);
A post processing can be applied to enhance the result of recovery. One can for instance average the result of NLmeans with the original, noisy, image.
% weighting factor lambda = .05; % averaging Mnl1 = (1-lambda)*Mnl0 + lambda*M;
Display.
clf; imageplot(clamp(Mnl0), ['NLmeans, SNR=' num2str(snr(M0,Mnl0)) 'dB'], 1,2,1); imageplot(clamp(Mnl1), ['Post-processing, SNR=' num2str(snr(M0,Mnl1)) 'dB'], 1,2,2);
Exercice 3: (the solution is exo3.m) Find the optimal weight lambda together with the optimal result Mnl1.
exo3;
Exercice 4: (the solution is exo4.m) Compare the non-local means denoising result with denoising using translation invariant wavelet hard thresholding.
exo4;
Exercice 5: (the solution is exo5.m) Explore the influence of the q and w parameters.
exo5;