Implemented BD3M

%Demo

clc

clearvars

close all

img = imread('Picture1.png');

%% BD3M Filter

sigma = 0.1225;

window_size = 8;

search_width = 19;

l2 = 0;

selection_number = 8;

l3 = 2.7;

img = imresize(img, [50 50]);

parfor channel = 1:3

noisy_image = padarray(img(:,:,channel),[search_width search_width], 'symmetric','both');

basic_result(:,:,channel) = first_step(noisy_image, sigma, window_size, search_width, l2, l3, selection_number);

basic_padded = padarray(basic_result(:,:,channel), [search_width search_width],'symmetric','both');

final_result(:,:,channel) = second_step(noisy_image,basic_padded, sigma, window_size, search_width, l2, selection_number);

end

final_result = uint8(final_result);

imshow(final_result)

function res = dct3(data, mode)

if(nargin > 2 && mode == 'inverse')

res = idct(data,'dim',3);

for i = 1:size(data,3)

res(:,:,i)=idct2(res(:,:,i));

end

else

res=zeros(size(data));

for i=1:size(data,2)

res(:,:,i)=dct2(data(:,:,i));

end

dct(res,'dim',3);

end

function basic_result = first_step(noisy_image, sigma, ws, sw, l2, l3, sn)

image_size = size(noisy_image);

numerator = zeros(size(noisy_image));

denominator = zeros(size(noisy_image));

bpr = (sw*2+1) - ws + 1; % number of blocks per row/col

center_block = bpr^2 / 2 + bpr/2 + 1;

for i = sw+1:image_size(1)-sw

for j = sw+1:image_size(2)-sw

window = noisy_image(i-sw:i+sw , j-sw:j+sw);

blocks = double(im2col(window, [ws ws], 'sliding'));

dist = zeros(size(blocks,2),1);

for k = 1:size(blocks,2)

tmp = wthresh(blocks(:,center_block),'h',sigma*l2)- ...

wthresh(blocks(:,k),'h',sigma*l2);

tmp = reshape(tmp, [ws ws]);

tmp = norm(tmp,2)^2;

dist(k) = tmp/(ws^2);

end

[~, inds] = sort(dist);

inds = inds(1:sn);

blocks = blocks(:, inds);

p = zeros([ws ws sn]);

for k = 1:sn

p(:,:,k) = reshape(blocks(:,k), [ws ws]);

end

p = dct3(p);

p = wthresh(p,'h',sigma*l3);

wp = 1/sum(p(:)>0);

p = dct3(p,'inverse');

for k = 1:sn

x = max(1,i-sw) + floor((center_block-1)/bpr);

y = max(1,j-sw) + (mod(center_block-1,bpr));

numerator(x:x+ws-1 , y:y+ws-1) = ...

numerator(x:x+ws-1 , y:y+ws-1) + (wp * p(:,:,k));

denominator(x:x+ws-1 , y:y+ws-1) = ...

denominator(x:x+ws-1 , y:y+ws-1) + wp;

end

basic_result = numerator./denominator;

basic_result = basic_result(sw+1:end-sw,sw+1:end-sw);

end

function result = second_step(noisy_image, basic_res, sigma, ws,sw, l2, sn)

image_size = size(noisy_image);

numerator = zeros(size(noisy_image));

denomerator = zeros(size(noisy_image));

bpr = (sw*2+1) - ws + 1; % number of blocks per row/col

center_block = bpr^2 / 2 + bpr/2 + 1;

for i = sw+1:image_size(1)-sw

for j = sw+1:image_size(2)-sw

window = noisy_image (i-sw:i+sw, j-sw:j+sw);

window2 = basic_res(i-sw:i+sw, j-sw:j+sw);

blocks = double(im2col(window, [ws ws], 'sliding'));

blocks2 = double(im2col(window2, [ws ws], 'sliding'));

dist = zeros(size(blocks,2),1);

for k = 1:size(blocks,2)

tmp = wthresh(blocks2(:,center_block),'h',sigma*l2)- ...

wthresh(blocks2(:,k),'h',sigma*l2);

tmp = reshape(tmp, [ws ws]);

tmp = norm(tmp,2)^2;

dist(k) = tmp/(ws^2);

end

[~, I] = sort(dist);

I = I(1:sn);

blocks = blocks(:, I);

blocks2 = blocks2(:, I);

p = zeros([ws ws sn]);

basic_p = zeros([ws ws sn]);

for k = 1:sn

p(:,:,k) = reshape(blocks(:,k), [ws ws]);

basic_p(:,:,k) = reshape(blocks2(:,k), [ws ws]);

end

basic_p = dct3(basic_p);

wp = zeros(sn,1);

for k = 1:sn

tmp = basic_p(:,:,k);

tmp = norm(tmp,1).^2;

wp(k) = tmp/(tmp+(sigma^2));

end

p = dct3(p);

for k = 1:sn

p(:,:,k) = p(:,:,k)*wp(k);

end

p = dct3(p,'inverse');

wp = 1/sum(wp(:).^2);

for k = 1:sn

q = p(:,:,k);

x = max(1,i-sw) + floor((center_block-1)/bpr);

y = max(1,j-sw) + (mod(center_block-1,bpr));

numerator(x:x+ws-1 , y:y+ws-1) = ...

numerator(x:x+ws-1 , y:y+ws-1) + (wp * q);

denomerator(x:x+ws-1 , y:y+ws-1) = ...

wp + denomerator(x:x+ws-1 , y:y+ws-1);

end

result = numerator./denomerator;

result = result(sw+1:end-sw,sw+1:end-sw);

end