Skip to content

Commit

Permalink
Sample code for I. Frosio, J. Kautz, Statistical Nearest Neighbors fo…
Browse files Browse the repository at this point in the history 
…r Image Denoising, IEEE Trans. Image Processing, 2018.
  • Loading branch information
@ifrosio
ifrosio committed Aug 31, 2018
1 parent 5ebb55c commit 3a18796
Show file tree
Hide file tree
Showing 19 changed files with 791 additions and 2 deletions.
65 changes: 65 additions & 0 deletions BilateralFilter/bilateralFilter.m
View file
@@ -0,0 +1,65 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function:img_f = bilateralFilter(img_n, halfFilterSize, sigma_r, sigma_s,
% sigma, offset)
% input: img_n, double image (single or multi channel) corrupted by iid
% Gaussian noise with std sigma;
% halfFilterSize, half size of the filter;
% sigma_s, sigma (spatial);
% sigma_r, sigma (range);
% sigma, noise std;
% offset, offset parameter (0 for NN, 1 for SNN).
% output: img_f, filtered image.
% scope: this function shows bilteral filter denoising with the NN / SNN
% sampling strategies; it is not designed to be computationally
% efficient.
% author: Iuri Frosio, ifrosio@nvidia.com
% ref: I. Frosio, J. Kautz, Statistical Neareast Neighbors for Image
% Denoising, IEEE Trans. Image Processing, 2018.
% license: Copyright (C) 2018 NVIDIA Corporation. All rights reserved.
% Licensed under the CC BY-NC-SA 4.0 license
% (https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function img_f = bilateralFilter(img_n, halfFilterSize, sigma_r, sigma_s, sigma, offset)

% init
filterSize = 2 * halfFilterSize + 1;
padded_img_n = padarray(img_n, [halfFilterSize halfFilterSize 0],'symmetric');
[ys, xs, cs] = size(img_n);
expected_squared_distance = 2 * sigma^2;
P = cs;
num = zeros(ys + 2 * halfFilterSize, xs + 2 * halfFilterSize, cs);
den = zeros(ys + 2 * halfFilterSize, xs + 2 * halfFilterSize, cs);
img_f = zeros(ys, xs, cs);

% spatial weights - these are fixed
[xxs,yys] = meshgrid(-halfFilterSize:halfFilterSize);
w_s = exp(-0.5 * (xxs.^2 + yys.^2) / (sigma_s^2));

% loop over all the pixels of the image
for y = halfFilterSize + 1 : ys + halfFilterSize
for x = halfFilterSize + 1 : ys + halfFilterSize

% compute sqaured distance from the central pixels
d2 = zeros(filterSize);
for c = 1 : cs
d2 = d2 + (padded_img_n(y-halfFilterSize:y+halfFilterSize, x-halfFilterSize:x+halfFilterSize, c) - ...
padded_img_n(y, x, c)).^2;
end
d2 = d2 / P;

% compute the range weights
w_r = exp(-0.5 * (abs(d2 - offset * expected_squared_distance)) / (sigma_r^2));

% weights
w = w_r .* w_s;
w = w / sum(w(:));

% filter
for c = 1 : cs
img_f(y - halfFilterSize, x - halfFilterSize, c) = sum(sum(w.*padded_img_n(y-halfFilterSize:y+halfFilterSize, x-halfFilterSize:x+halfFilterSize, c)));
end
end
end


end
61 changes: 61 additions & 0 deletions BilateralFilter/bilateralFilterTest.m
View file
@@ -0,0 +1,61 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: bilateralTest
% input: none
% output: none
% scope: this function shows how to call bilateralFilter.m. It gives an
% example of denoising using the bilateral filter with the
% NN / SNN sampling strategies.
% author: Iuri Frosio, ifrosio@nvidia.com
% ref: I. Frosio, J. Kautz, Statistical Neareast Neighbors for Image
% Denoising, IEEE Trans. Image Processing, 2018.
% license: Copyright (C) 2018 NVIDIA Corporation. All rights reserved.
% Licensed under the CC BY-NC-SA 4.0 license
% (https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function bilateralTest

% parameters
sigma = 0.3; % noise std
halfFilterSize = 5; % half size of the patch
sigma_s = 3; % spatial sigma
sigma_r = sigma; % range sigma (assumed equal to the noise std)

% create an image to denoise
[x,y] = meshgrid(1:512);
img = sin(1.4*x/512*pi) + cos((1.3*(x+y)/512*pi).^2) - cos((0.34*(2*x+y)/512*pi).^4);
img = (img - min(img(:)));
img = img/max(img(:));
img = imresize(img, 0.25);
img(:,:,2) = 0.25 + 0.25*img(:,:,1);
img(:,:,3) = 0.9 - 0.7* img(:,:,1);
img_n = img + randn(size(img)) * sigma;

% denoise (nn, offset = 0)
offset = 0;
img_f_nn = bilateralFilter(img_n, halfFilterSize, sigma_r, sigma_s, sigma, offset);

% denoise (snn, offset = 0.8)
offset = 0.8;
img_f_snn = bilateralFilter(img_n, halfFilterSize, sigma_r, sigma_s, sigma, offset);

% errors
mse_n = mean((img(:)-img_n(:)).^2);
mse_nn = mean((img(:)-img_f_nn(:)).^2);
mse_snn = mean((img(:)-img_f_snn(:)).^2);

figure(1);
clf;
subplot(221);
imshow(img);
title('Image');
subplot(222);
imshow(img_n);
title(['Noisy image - MSE = ' num2str(mse_n)]);
subplot(223);
imshow(img_f_nn);
title(['Filtered image [NN] - MSE = ' num2str(mse_nn)]);
subplot(224);
imshow(img_f_snn);
title(['Filtered image [SNN] - MSE = ' num2str(mse_snn)]);

end
4 changes: 4 additions & 0 deletions BilateralFilter/runme.m
View file
@@ -0,0 +1,4 @@
% run me!
function runme
bilateralFilterTest;
end
110 changes: 110 additions & 0 deletions NLM/nlm.m
View file
@@ -0,0 +1,110 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: img_f = nlm(img_n, halfPatchSize, windowHalfSearchSize, N_n,
% sigma, h, offset)
% input: img_n, double image (single or multi channel) corrupted by iid
% Gaussian noise with std sigma;
% halfPatchSize, half size of the patch used by NLM;
% windowHalfSearchSize, half size of the windows used to search
% for neighbors;
% N_n, number of neighbors used by NLM;
% sigma, std of the noise;
% h, filtering parameter for NLM;
% offset, offset parameter (0 for NN, 1 for SNN).
% output: img_f, filtered image.
% scope: this function shows NLM denoising with the NN / SNN sampling
% strategies; it is not designed to be computationally efficient.
% author: Iuri Frosio, ifrosio@nvidia.com
% ref: I. Frosio, J. Kautz, Statistical Neareast Neighbors for Image
% Denoising, IEEE Trans. Image Processing, 2018.
% license: Copyright (C) 2018 NVIDIA Corporation. All rights reserved.
% Licensed under the CC BY-NC-SA 4.0 license
% (https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function img_f = nlm(img_n, halfPatchSize, windowHalfSearchSize, N_n, sigma, h, offset)

% init
[xs, ys, cs] = size(img_n);
patchSize = 2 * halfPatchSize + 1;
P = cs * patchSize^2;
expected_squared_distance = 2 * sigma^2;

% Init buffers
neighbors_indexes = zeros(ys, xs, N_n);
neighbors_d2 = ones(ys, xs, N_n) * inf;
padded_img_n = padarray(img_n, [windowHalfSearchSize windowHalfSearchSize 0], 'symmetric');

% For each shift
index = 0;
for dy = -windowHalfSearchSize:windowHalfSearchSize
for dx = -windowHalfSearchSize:windowHalfSearchSize

% shift index
index = index + 1;

% shifted image (with mirroring)
shifted_img_n = padded_img_n((windowHalfSearchSize + 1 + dy):(windowHalfSearchSize + ys + dy), ...
(windowHalfSearchSize + 1 + dx):(windowHalfSearchSize + xs + dx), ...
:);

% squared distance between the shifted and the reference image
current_d2 = imfilter(sum((shifted_img_n - img_n).^2, 3), ones(patchSize)/P, 'symmetric');
current_indexes = ones(ys, xs) * index;

% update neighbors
for n = 1 : N_n

% is the current neighbor closer than the stored one?
swap = abs(current_d2 - offset * expected_squared_distance) < abs(neighbors_d2(:, :, n) - offset * expected_squared_distance);

% swap (indexes, distances)
neighbors_indexes_n = neighbors_indexes(:, :, n);
buffer_indexes_n = neighbors_indexes_n;
neighbors_indexes_n(swap) = current_indexes(swap);
current_indexes(swap) = buffer_indexes_n(swap);
neighbors_indexes(:, :, n) = neighbors_indexes_n;

neighbors_d2_n = neighbors_d2(:, :, n);
buffer_d2_n = neighbors_d2_n;
neighbors_d2_n(swap) = current_d2(swap);
current_d2(swap) = buffer_d2_n(swap);
neighbors_d2(:, :, n) = neighbors_d2_n;

end
end
end

% init num / den
num = zeros(ys, xs, cs);
den = zeros(ys, xs, cs);

% do another loop on the possible neighbors to filter
index = 0;
for dy = -windowHalfSearchSize:windowHalfSearchSize
for dx = -windowHalfSearchSize:windowHalfSearchSize

% shift index
index = index + 1;

% shifted image (with mirroring)
shifted_img_n = padded_img_n((windowHalfSearchSize + 1 + dy):(windowHalfSearchSize + ys + dy), ...
(windowHalfSearchSize + 1 + dx):(windowHalfSearchSize + xs + dx), ...
:);

% For every neighbors
for n = 1 : N_n

% weights
buffer_weights = (neighbors_indexes(:, :, n) == index) .* exp(-(max(0, neighbors_d2(:, :, n) - 2 * sigma^2))/(h^2));
weights = repmat(imfilter(buffer_weights, ones(patchSize), 'symmetric'), [1 1 cs]);
num = num + weights .* shifted_img_n;
den = den + weights;

end
end
end

% filtered image
img_f = num./den;

end
63 changes: 63 additions & 0 deletions NLM/nlmTest.m
View file
@@ -0,0 +1,63 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: nlmTest
% input: none
% output: none
% scope: this function shows how to call nlm.m. It gives an example of
% denoising using the NLM filter with the NN / SNN sampling
% strategies.
% author: Iuri Frosio, ifrosio@nvidia.com
% ref: I. Frosio, J. Kautz, Statistical Neareast Neighbors for Image
% Denoising, IEEE Trans. Image Processing, 2018.
% license: Copyright (C) 2018 NVIDIA Corporation. All rights reserved.
% Licensed under the CC BY-NC-SA 4.0 license
% (https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nlmTest

% parameters
sigma = 0.3; % noise std
halfPatchSize = 3; % half size of the patch
windowSearchHalfSize = 6; % half size for searching the neighbors
N_n = 16; % numer of neighbors
h = 0.3 * sigma^2; % nlm filtering parameter

% create an image to denoise
[x,y] = meshgrid(1:512);
img = sin(1.4*x/512*pi) + cos((1.3*(x+y)/512*pi).^2) - cos((0.34*(2*x+y)/512*pi).^4);
img = (img - min(img(:)));
img = img/max(img(:));
img = imresize(img, 0.25);
img(:,:,2) = 0.25 + 0.25*img(:,:,1);
img(:,:,3) = 0.9 - 0.7* img(:,:,1);
img_n = img + randn(size(img)) * sigma;


% denoise (nn, offset = 0)
offset = 0;
img_f_nn = nlm(img_n, halfPatchSize, windowSearchHalfSize, N_n, sigma, h, offset);

% denoise (snn, offset = 0.8)
offset = 0.8;
img_f_snn = nlm(img_n, halfPatchSize, windowSearchHalfSize, N_n, sigma, h, offset);

% errors
mse_n = mean((img(:)-img_n(:)).^2);
mse_nn = mean((img(:)-img_f_nn(:)).^2);
mse_snn = mean((img(:)-img_f_snn(:)).^2);

figure(1);
clf;
subplot(221);
imshow(img);
title('Image');
subplot(222);
imshow(img_n);
title(['Noisy image - MSE = ' num2str(mse_n)]);
subplot(223);
imshow(img_f_nn);
title(['Filtered image [NN] - MSE = ' num2str(mse_nn)]);
subplot(224);
imshow(img_f_snn);
title(['Filtered image [SNN] - MSE = ' num2str(mse_snn)]);

end
4 changes: 4 additions & 0 deletions NLM/runme.m
View file
@@ -0,0 +1,4 @@
% run me!
function runme
nlmTest;
end
23 changes: 21 additions & 2 deletions README.md
View file
@@ -1,2 +1,21 @@
# SNN
Matlab code implementation the modified Non Local Means and Bilateral filters, as described in I. Frosio, J. Kautz, Statistical Nearest Neighbors for Image Denoising, IEEE Trans. Image Processing, 2018. The repository also includes the Matlab code to replicate the results of the toy problem described in the paper.
# SNN (Statistical Nearest Neighbors)

Matlab code implementation of the modified Non Local Means and Bilateral filters, as described in I. Frosio, J. Kautz, Statistical Nearest Neighbors for Image Denoising, IEEE Trans. Image Processing, 2018. The repository also includes the Matlab code to replicate the results of the toy problem described in the paper.

# Usage

**/NLM/nlm.m**: Matlab code to filter a user defined image thorugh the modified Non Local Means filter (accordinagly to the SNN approach) and user defined filtering parameters. The traditional Non Local Means filter is achieved passing offset = 0 in input.

**/NLM/runme.m**: example usage of /NLM/nlm.m.

**/BilateralFilter/bilateralFilter.m**: Matlab code to filter a user defined image thorugh the modified Bilateral filter (accordinagly to the SNN approach) and user defined filtering parameters. The traditional Bilateral filter is achieved passing offset = 0 in input.

**/BilateralFilter/runme.m**: example usage of /BilateralFilter/bilateralFilter.m.

**/ToyProblem/toyProblem.m**: Matlab code to replicate the results on the toy problem described in our paper.

**/ToyProblem/runme.m**: example usage of ToyProblem/toyproblem.m.

# References

I. Frosio, J. Kautz, Statistical Nearest Neighbors for Image Denoising, IEEE Trans. Image Processing, 2018.
15 changes: 15 additions & 0 deletions ToyProblem/Phi_.m
View file
@@ -0,0 +1,15 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: Phix = Phi(x)
% input:
% output:
% scope:
% author: Iuri Frosio, ifrosio@nvidia.com
% ref: I. Frosio, J. Kautz, Statistical Neareast Neighbors for Image
% Denoising, IEEE Trans. Image Processing, 2018.
% license: Copyright (C) 2018 NVIDIA Corporation. All rights reserved.
% Licensed under the CC BY-NC-SA 4.0 license
% (https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Phix = Phi(x)

Phix = 0.5 * (1 + erf(x/sqrt(2)));
15 changes: 15 additions & 0 deletions ToyProblem/Phi_mix_.m
View file
@@ -0,0 +1,15 @@
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function: y = Phi_mix_(x, p_r, mu, sigma, p_f, mu_f, sigma_f)
% input:
% output:
% scope:
% author: Iuri Frosio, ifrosio@nvidia.com
% ref: I. Frosio, J. Kautz, Statistical Neareast Neighbors for Image
% Denoising, IEEE Trans. Image Processing, 2018.
% license: Copyright (C) 2018 NVIDIA Corporation. All rights reserved.
% Licensed under the CC BY-NC-SA 4.0 license
% (https://creativecommons.org/licenses/by-nc-sa/4.0/legalcode).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = Phi_mix_(x, p_r, mu, sigma, p_f, mu_f, sigma_f)
y = p_r * Phi_((x - mu)/sigma) + p_f * Phi_((x - mu_f) / sigma_f);
end

0 comments on commit 3a18796

Please sign in to comment.