BM3D denoise

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Copyright � 2005 Tampere University of Technology. All rights reserved. % This work should only be used for nonprofit purposes. % % AUTHORS: % Kostadin Dabov (2005), email: kostadin.dabov _at_ tut.fi % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % 1. OVERVIEW % This script implements the BM3DDFT algorithm for attenuation of % additive white Gaussian noise from images. BM3DDFT stands for "image % denoising algorithm with block-matching and filtering in 3D DFT % domain". % % 2. USAGE % 1) Set the filename of a noise-free image (parameter "image_name") % 2) Set the desired standard deviation of the AWGN (parameter "sigma") % 3) Invoke the script from MATLAB's command prompt % % 3. REQUIREMENTS % The compiled files work on Linux or Windows, Pentium-based machines. % The script is tested to work with Matlab v.7. % % 4. REFERENCE % [1] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image denoising % with block-matching and 3D filtering," in Electronic Imaging�06, % Proc. SPIE 6064, no. 6064A-30, San Jose, California USA, 2006. % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Select an image filename (might contain path also) image_name = [ % 'boat.png' % 'lena.png' 'house.png' % 'barbara.png' % 'peppers256.png' % 'couple.tif' % 'hill.tif' % 'man.tif' ]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Standard deviation of the white Gaussian noise (WGN) sigma = 25; %% standard deviation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Below are the parameters that reproduce the results reported in [1]. %%%% Further information for the parameters can be found the same paper. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Hard-thresholding (gives initial estimate) parameters Nstep = 4; %% sliding step to process every next refernece block N2 = 28; %% maximum number of similar blocks (maximum 3rd dimension size of a 3D block) Ns = 73; %% search window's length of the side beta = 4; %% parameter of the 2D Kaiser window tau_match = 0.233; %% threshold for the block distance (d-distance) lambda_thr2D = 0.82; %% threshold parameter for the coarse initial denoising used in the d-distance measure lambda_thr3D = 0.75; %% threshold parameter for the hard-thresholding in 3D DFT domain %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Wiener filtering (gives final estimate) parameters Nstep_wiener = 3; N2_wiener = 72; Ns_wiener = 35; beta_wiener = 3; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Parameters selected automatically based on sigma N1 = min(max(round(sigma*0.08 + 9.4), 8), 13); %% block-size used for hard-thresholding N1 = N1 - 1 + rem(N1,2); %% ensure block size of odd size N1_wiener = min(max(round(0.15*sigma + 5.625), 7), 11); %% block-size used for wiener filtering N1_wiener = N1_wiener - 1 + rem(N1_wiener,2); %% ensure block size of odd size tau_match_wiener = sigma/4000 + 0.0105; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Note: touch below this point only if you know what you are doing! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Make parameters compatible with the interface of the mex-functions sigma = sigma/255; % put in range [0,1] Ns_half = (Ns-1)/2; % half the search window size needed for the hard-thresholding routine Ns_wiener_half = (Ns_wiener-1)/2; % half the search window size needed for the wiener routine tau_match = tau_match*N1; % because N1 is constant in the d-discatnce formula, % we pre-multiply tau_match with it. tau_match_wiener = tau_match_wiener*N1_wiener; % because N1_wiener is constant in the % discatnce formula for Wiener block-matching, % we pre-multiply tau_match_wiener with it. Wwin2D = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the hard-thresholding part Wwin2D_wiener = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the Wiener filtering part %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Read an image, generate noise and add it to the image y = im2double(imread(image_name)); %% read a noise-free image randn('seed', 0); %% generate seed z = y + sigma*randn(size(y)); %% create a noisy image [Xv, Xh] = size(y); %% obtain image sizes %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Print image information to the screen fprintf(sprintf('Image: %s (%dx%d), sigma: %.1f\n', image_name, Xv, Xh, sigma*255)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Initial estimate by hard-thresholding filtering tic; AVG = denfft2d3d_c(z, 'unused_arg', Nstep, N1, N2, lambda_thr2D,... lambda_thr3D, tau_match, Ns_half, sigma, Wwin2D); estimate_elapsed_time = toc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Calculate initial estimate's PSNR and ISNR, and print them PSNR_NOISY = 10*log10(1/mean2((y-z).^2)); PSNR_INITIAL_ESTIMATE = 10*log10(1/mean2((y-AVG).^2)); ISNR_INITIAL_ESTIMATE = PSNR_INITIAL_ESTIMATE - PSNR_NOISY; fprintf(sprintf('INITIAL ESTIMATE, ISNR: %.2f dB, PSNR: %.2f dB\n', ... ISNR_INITIAL_ESTIMATE, PSNR_INITIAL_ESTIMATE)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Final estimate by Wiener filtering (using the hard-thresholding initial estimate) tic; AVGW = denfft2d3d_wiener_c(z, 'unused_arg', Nstep_wiener, N1_wiener, N2_wiener, 'unused_arg', ... 'unused_arg', tau_match_wiener, Ns_wiener_half, sigma, Wwin2D_wiener, AVG); wiener_elapsed_time = toc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Calculate final estimate's PSNR and ISNR, print them, and show the %%%% denoised image next to the noisy one PSNR_FINAL_ESTIMATE = 10*log10(1/mean2((y-AVGW).^2)); ISNR_FINAL_ESTIMATE = PSNR_FINAL_ESTIMATE - PSNR_NOISY; fprintf(sprintf('FINAL ESTIMATE (total time: %.1f sec), ISNR: %.2f dB, PSNR: %.2f dB\n', ... wiener_elapsed_time + estimate_elapsed_time, ISNR_FINAL_ESTIMATE, PSNR_FINAL_ESTIMATE)); imshow([z, AVGW]); title(sprintf('Noisy (standard deviation %d) and denoised (PSNR %.2f dB) images for "%s"', sigma*255, PSNR_FINAL_ESTIMATE, image_name)); return;