Non Local Means

This numerical tour study image denoising using non-local means.

Installing toolboxes and setting up the path.
Patches in Images
Dimensionality Reduction with PCA
Non-local Filter
Localizing the Non-local Means
Non-local Means for Denoising

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 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/');

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;