BM3D/VBM3D.m

   1 function [PSNR_FINAL_ESTIMATE, y_hat_wi] = VBM3D(Xnoisy, sigma, NumberOfFrames, dump_information, Xorig, bm3dProfile)
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %
   4 %  VBM3D is a Matlab function for attenuation of additive white Gaussian
   5 %  noise from grayscale videos. This algorithm reproduces the results from the article:
   6 %
   7 %  [1] K. Dabov, A. Foi, and K. Egiazarian, "Video denoising by sparse 3D
   8 %  transform-domain collaborative filtering," European Signal Processing
   9 %  Conference (EUSIPCO-2007), September 2007. (accepted)
  10 %
  11 %  INTERFACE:
  12 %
  13 %  [PSNR, Xest] = VBM3D(Xnoisy, Sigma, NFrames, PrintInfo, Xorig)
  14 %
  15 %  INPUTS:
  16 %   1)  Xnoisy     --> A filename of a noisy .avi video, e.g. Xnoisy = 'gstennisg20.avi'
  17 %        OR
  18 %       Xnoisy     --> A 3D matrix of a noisy video in a  (floating point data in range [0,1],
  19 %                                                     or in [0,255])
  20 %   2)  Sigma --> Noise standard deviation (assumed range is [0,255], no matter what is
  21 %                                           the input's range)
  22 %
  23 %   3)  NFrames   (optional paremter!) --> Number of frames to process. If set to 0 or
  24 %                                          ommited, then process all frames (default: 0).
  25 %
  26 %   4)  PrintInfo (optional paremter!) --> If non-zero, then print to screen and save
  27 %                                          the denoised video in .AVI
  28 %                                          format. (default: 1)
  29 %
  30 %   5)  Xorig     (optional paremter!) --> Original video's filename or 3D matrix
  31 %                                          If provided, PSNR, ISNR will be computed.
  32 %
  33 %   NOTE: If Xorig == Xnoisy, then artificial noise is added internally and the
  34 %   obtained noisy video is denoised.
  35 %
  36 %  OUTPUTS:
  37 %
  38 %   1) PSNR --> If Xorig is valid video, then this contains the PSNR of the
  39 %                denoised one
  40 %
  41 %   1) Xest --> Final video estimate in a 3D matrix (intensities in range [0,1])
  42 %
  43 %   *) If "PrintInfo" is non-zero, then save the denoised video in the current
  44 %       MATLAB folder.
  45 %
  46 %  USAGE EXAMPLES:
  47 %
  48 %     1) Denoise a noisy (clipped in [0,255] range) video sequence, e.g.
  49 %        'gsalesmang20.avi' corrupted with AWGN with std. dev. 20:
  50 %
  51 %          Xest = VBM3D('gsalesmang20.avi', 20, 0, 1);
  52 %
  53 %     2) The same, but also print PSNR, ISNR numbers.
  54 %
  55 %          Xest = VBM3D('gsalesmang20.avi', 20, 0, 1, 'gsalesman.avi');
  56 %
  57 %     3) Add artificial noise to a video, then denoise it (without
  58 %        considering clipping in [0,255]):
  59 %
  60 %          Xest = VBM3D('gsalesman.avi', 20, 0, 1, 'gsalesman.avi');
  61 %
  62 %
  63 %  RESTRICTIONS:
  64 %
  65 %     Since the video sequences are read into memory as 3D matrices,
  66 %     there apply restrictions on the input video size, which are thus
  67 %     proportional to the maximum memory allocatable by Matlab.
  68 %
  69 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  70 %
  71 % Copyright © 2007 Tampere University of Technology. All rights reserved.
  72 % This work should only be used for nonprofit purposes.
  73 %
  74 % AUTHORS:
  75 %     Kostadin Dabov, email: dabov _at_ cs.tut.fi
  76 %
  77 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  78
  79 % If no input argument is provided, then use these internal ones:
  80 if exist('sigma', 'var') ~= 1,
  81     Xnoisy = 'gsalesmang20.avi'; Xorig = 'gsalesman.avi'; sigma = 20;
  82     %Xnoisy = 'gstennisg20.avi';  Xorig = 'gstennis.avi';  sigma = 20;
  83     %Xnoisy = 'gflowersg20.avi';   Xorig = 'gflower.avi';   sigma = 20;
  84
  85     %Xnoisy = 'gsalesman.avi'; Xorig = Xnoisy; sigma = 20;
  86
  87     NumberOfFrames = 0; %% 0 means process ALL frames.
  88 end
  89
  90
  91
  92 if exist('dump_information', 'var') ~= 1,
  93     dump_information = 1; % 1 -> print informaion to the screen and save the processed video as an AVI file
  94 end
  95
  96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  97 %%%% Obtain infromation about the input noisy video
  98 %%%%
  99 if (ischar(Xnoisy) == 1), % if the input is a video filename
 100     isCharacterName = 1;
 101     Xnoisy_name = Xnoisy;
 102     videoInfo = aviinfo(Xnoisy);
 103     videoHeight = videoInfo.Height;
 104     videoWidth = videoInfo.Width;
 105     TotalFrames = videoInfo.NumFrames;
 106 elseif length(size(Xnoisy)) == 3% the input argument is a 3D video (spatio-temporal) matrix
 107     Xnoisy_name = 'Input 3D matrix';
 108     isCharacterName = 0;
 109     [videoHeight, videoWidth, TotalFrames] = size(Xnoisy);
 110 else
 111     fprintf('Oops! The input argument Xnoisy should be either a filename or a 3D matrix!\n');
 112     PSNR_FINAL_ESTIMATE = 0;
 113     y_hat_wi = 0;
 114     return;
 115 end
 116
 117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 118 %%%% Check if we want to process all frames, and save as 'NumberOfFrames'
 119 %%%% the desired number of frames to process
 120 %%%%
 121 if exist('NumberOfFrames', 'var') == 1,
 122     if NumberOfFrames <= 0,
 123         NumberOfFrames = TotalFrames;
 124     else
 125         NumberOfFrames = max(min(NumberOfFrames, TotalFrames), 1);
 126     end
 127 else
 128     NumberOfFrames = TotalFrames;
 129 end
 130
 131 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 132 %%%%  Quality/complexity trade-off
 133 %%%%
 134 %%%%  'np' --> Normal Profile (balanced quality)
 135 %%%%  'lc' --> Low Complexity Profile (fast, lower quality)
 136 %%%%
 137 if (exist('bm3dProfile', 'var') ~= 1)
 138     bm3dProfile         = 'np';
 139 end
 140
 141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 142 %%%% Parameters for the Normal Profile.
 143 %%%%
 144 %%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
 145 transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
 146 transform_2D_Wiener_name = 'dct';     %% transform used for the Wiener filt. of size N1_wiener x N1_wiener
 147 transform_3rd_dim_name   = 'haar'; %% tranform used in the 3-rd dim, the same for HT and Wiener filt.
 148
 149 %%%% Step 1: Hard-thresholding (HT) parameters:
 150 denoiseFrames       = min(9, NumberOfFrames); % number of frames in the temporalwindow (should not exceed the total number of frames 'NumberOfFrames')
 151 N1                  = 8;  %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
 152 Nstep               = 6;  %% sliding step to process every next refernece block
 153 N2                  = 8;  %% maximum number of similar blocks (maximum size of the 3rd dimension of the 3D groups)
 154 Ns                  = 7;  %% length of the side of the search neighborhood for full-search block-matching (BM)
 155 Npr                 = 5;  %% length of the side of the motion-adaptive search neighborhood, use din the predictive-search BM
 156 tau_match           = 3000; %% threshold for the block distance (d-distance)
 157 lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D DFT domain
 158 dsub                = 7;  %% a small value subtracted from the distnce of blocks with the same spatial coordinate as the reference one
 159 Nb                  = 2;  %% number of blocks to follow in each next frame, used in the predictive-search BM
 160 beta                = 2.0; %% the beta parameter of the 2D Kaiser window used in the reconstruction
 161
 162 %%%% Step 2: Wiener filtering parameters:
 163 denoiseFramesW      = min(9, NumberOfFrames);
 164 N1_wiener           = 7;
 165 Nstep_wiener        = 4;
 166 N2_wiener           = 8;
 167 Ns_wiener           = 7;
 168 Npr_wiener          = 5;
 169 tau_match_wiener    = 1500;
 170 beta_wiener         = 2.0;
 171 dsub_wiener         = 3;
 172 Nb_wiener           = 2;
 173
 174 %%%% Block-matching parameters:
 175 stepFS              = 1; %% step that forces to switch to full-search BM, "1" implies always full-search
 176 smallLN             = 3; %% if stepFS > 1, then this specifies the size of the small local search neighb.
 177 stepFSW             = 1;
 178 smallLNW            = 3;
 179 thrToIncStep        = 8;  %% used in the HT filtering to increase the sliding step in uniform regions
 180
 181
 182 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 183 %%%% Parameters for the Low Complexity Profile.
 184 %%%%
 185 if strcmp(bm3dProfile, 'lc') == 1,
 186     lambda_thr3D = 2.8;
 187     smallLN   = 2;
 188     smallLNW  = 2;
 189     denoiseFrames  = min(5, NumberOfFrames);
 190     denoiseFramesW = min(5, NumberOfFrames);
 191     N2_wiener = 4;
 192     N2 = 4;
 193     Ns = 3;
 194     Ns_wiener = 3;
 195     NB = 1;
 196     Nb_wiener = 1;
 197 end
 198
 199 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 200 %%%% Parameters for the High Profile.
 201 %%%%
 202 if strcmp(bm3dProfile, 'hi') == 1,
 203     Nstep        = 3;
 204     Nstep_wiener = 3;
 205 end
 206
 207 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 208 %%%% Parameters for the "Very Noisy" Profile.
 209 %%%%
 210 if sigma > 30,
 211     N1 = 8;
 212     N1_wiener = 8;
 213     Nstep = 6;
 214     tau_match    = 4500;
 215     tau_match_wiener    = 3000;
 216 end
 217
 218
 219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 221 %%%% Note: touch below this point only if you know what you are doing!
 222 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 223 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 224
 225 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 226 %%%% Extract the input noisy video and make sure intensities are in [0,1]
 227 %%%% interval, using single-precision float
 228 if isCharacterName,
 229     mno = aviread(Xnoisy_name);
 230     z = zeros([videoHeight, videoWidth, NumberOfFrames], 'single');
 231     for cf = 1:NumberOfFrames
 232         z(:,:,cf) = single(mno(cf).cdata(:,:,1)) * 0.0039216; % 1/255 = 0.0039216
 233     end
 234     clear  mno
 235 else
 236     if isinteger(Xnoisy) == 1,
 237         z = single(Xnoisy) * 0.0039216; % 1/255 = 0.0039216
 238     elseif isfloat(Xnoisy) == 0,
 239         fprintf('Unknown format of "Xnoisy"! Must be a filename (array of char) or a 3D array of either floating point data (range [0,1]) or integer data (range [0,255]). \n');
 240         return;
 241     else
 242         z = single(Xnoisy);
 243     end
 244 end
 245
 246 clear Xnoisy;
 247
 248 %%%% If the original video is provided, then extract it to 'Xorig'
 249 %%%% which is later used to compute PSNR and ISNR
 250 if exist('Xorig', 'var') == 1,
 251     randn('seed', 0);
 252     if ischar(Xorig) == 0,
 253         if isinteger(Xorig) == 1,
 254             y = single(Xorig) * 0.0039216; % 1/255 = 0.0039216
 255         elseif isfloat(Xorig) == 0,
 256             fprintf('Unknown format of "Xorig"! Must be a filename (array of char) or a 3D array of either floating point data (range [0,1]) or integer data (range [0,255]). \n');
 257             return;
 258         else
 259             y = single(Xorig);
 260         end
 261     else
 262         if strcmp(Xorig, Xnoisy_name) == 1, %% special case, noise is aritifically added
 263             y = z;
 264             z = z + (sigma/255) * randn(size(z));
 265         else
 266             mo = aviread(Xorig);
 267             y = zeros([videoHeight, videoWidth, NumberOfFrames], 'single');
 268             for cf = 1:NumberOfFrames
 269                 y(:,:,cf) = single(mo(cf).cdata(:,:,1)) * 0.0039216; % 1/255 = 0.0039216
 270             end
 271             clear mo
 272         end
 273     end
 274 end
 275
 276 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 277 %%%% Create transform matrices, etc.
 278 %%%%
 279 decLevel           = 0;    %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
 280 decLevel3          = 0;    %% dec. level for the wavelet transform in the 3rd dimension
 281
 282 [Tfor, Tinv]       = getTransfMatrix(N1, transform_2D_HT_name, decLevel); %% get (normalized) forward and inverse transform matrices
 283 [TforW, TinvW]     = getTransfMatrix(N1_wiener, transform_2D_Wiener_name); %% get (normalized) forward and inverse transform matrices
 284 thr_mask           = ones(N1); %% N1xN1 mask of threshold scaling coeff. --- by default there is no scaling, however the use of different thresholds for different wavelet decompoistion subbands can be done with this matrix
 285
 286 if (strcmp(transform_3rd_dim_name, 'haar') == 1 || strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
 287     %%% Fast internal transform is used, no need to generate transform
 288     %%% matrices.
 289     hadper_trans_single_den         = {};
 290     inverse_hadper_trans_single_den = {};
 291 else
 292     %%% Create transform matrices. The transforms are later computed by
 293     %%% matrix multiplication with them
 294     for hh = [1 2 4 8 16 32];
 295         [Tfor3rd, Tinv3rd]   = getTransfMatrix(hh, transform_3rd_dim_name, decLevel3);
 296         hadper_trans_single_den{hh}         = single(Tfor3rd);
 297         inverse_hadper_trans_single_den{hh} = single(Tinv3rd');
 298     end
 299 end
 300
 301 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 302 %%%% 2D Kaiser windows that scale the reconstructed blocks
 303 %%%%
 304 if beta_wiener==2 & beta==2 & N1_wiener==7 & N1==8 % hardcode the window function so that the signal processing toolbox is not needed by default
 305     Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
 306         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 307         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 308         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 309         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 310         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 311         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 312         0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924 ];
 313     Wwin2D_wiener = [ 0.1924    0.3151    0.4055    0.4387    0.4055    0.3151    0.1924;
 314         0.3151    0.5161    0.6640    0.7184    0.6640    0.5161    0.3151;
 315         0.4055    0.6640    0.8544    0.9243    0.8544    0.6640    0.4055;
 316         0.4387    0.7184    0.9243    1.0000    0.9243    0.7184    0.4387;
 317         0.4055    0.6640    0.8544    0.9243    0.8544    0.6640    0.4055;
 318         0.3151    0.5161    0.6640    0.7184    0.6640    0.5161    0.3151;
 319         0.1924    0.3151    0.4055    0.4387    0.4055    0.3151    0.1924 ];
 320 else
 321     Wwin2D           = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
 322     Wwin2D_wiener    = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the aggregation of the Wiener filt. part
 323 end
 324
 325 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 326 %%%% Read an image, generate noise and add it to the image
 327 %%%%
 328
 329 l2normLumChrom = ones(NumberOfFrames,1); %%% NumberOfFrames == nSl !
 330
 331 if dump_information == 1,
 332     fprintf('Video: %s (%dx%dx%d), sigma: %.1f\n', Xnoisy_name, videoHeight, videoWidth, NumberOfFrames, sigma);
 333 end
 334
 335 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 336 %%%% Initial estimate by hard-thresholding filtering
 337 tic;
 338 y_hat = bm3d_thr_video(z, hadper_trans_single_den, Nstep, N1, N2, 0,...
 339     lambda_thr3D, tau_match*N1*N1/(255*255), (Ns-1)/2, sigma/255, thrToIncStep, single(Tfor), single(Tinv)', inverse_hadper_trans_single_den, single(thr_mask), 'unused arg', dsub*dsub/255, l2normLumChrom, Wwin2D, (Npr-1)/2, stepFS, denoiseFrames, Nb );
 340 estimate_elapsed_time = toc;
 341
 342 if exist('Xorig', 'var') == 1,
 343     PSNR_INITIAL_ESTIMATE = 10*log10(1/mean((double(y(:))-double(y_hat(:))).^2));
 344     PSNR_NOISE = 10*log10(1/mean((double(y(:))-double(z(:))).^2));
 345     ISNR_INITIAL_ESTIMATE = PSNR_INITIAL_ESTIMATE - PSNR_NOISE;
 346
 347     if dump_information == 1,
 348         fprintf('BASIC ESTIMATE (time: %.1f sec), PSNR: %.3f dB, ISNR: %.3f dB\n', ...
 349             estimate_elapsed_time, PSNR_INITIAL_ESTIMATE, ISNR_INITIAL_ESTIMATE);
 350     end
 351 end
 352
 353
 354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 355 % %%%% Final estimate by Wiener filtering (using the hard-thresholding
 356 % initial estimate)
 357 tic;
 358 y_hat_wi = bm3d_wiener_video(z, y_hat, hadper_trans_single_den, Nstep_wiener, N1_wiener, N2_wiener, ...
 359     'unused_arg', tau_match_wiener*N1_wiener*N1_wiener/(255*255), (Ns_wiener-1)/2, sigma/255, 'unused arg', single(TforW), single(TinvW)', inverse_hadper_trans_single_den, 'unused arg', dsub_wiener*dsub_wiener/255, l2normLumChrom, Wwin2D_wiener, (Npr_wiener-1)/2, stepFSW, denoiseFramesW, Nb_wiener );
 360
 361 % In case the input noisy video is clipped in [0,1], then apply declipping
 362 if isCharacterName
 363     if exist('Xorig', 'var') == 1
 364         if ~strcmp(Xorig, Xnoisy_name)
 365             [y_hat_wi] = ClipComp16b(sigma/255, y_hat_wi);
 366         end
 367     else
 368         [y_hat_wi] = ClipComp16b(sigma/255, y_hat_wi);
 369     end
 370 end
 371
 372 wiener_elapsed_time = toc;
 373
 374
 375
 376 PSNR_FINAL_ESTIMATE = 0;
 377 if exist('Xorig', 'var') == 1,
 378     PSNR_FINAL_ESTIMATE = 10*log10(1/mean((double(y(:))-double(y_hat_wi(:))).^2));
 379     ISNR_FINAL_ESTIMATE = PSNR_FINAL_ESTIMATE - 10*log10(1/mean((double(y(:))-double(z(:))).^2));
 380 end
 381
 382 if dump_information == 1,
 383
 384     text_psnr = '';
 385     if exist('Xorig', 'var') == 1
 386
 387
 388         %%%% Un-comment the following to print the PSNR of each frame
 389         %
 390         %     PSNRs = zeros(NumberOfFrames,1);
 391         %     for ii = [1:NumberOfFrames],
 392         %         PSNRs(ii) = 10*log10(1/mean2((y(:,:,ii)-y_hat_wi(:,:,ii)).^2));
 393         %         fprintf(['Frame: ' sprintf('%d',ii) ', PSNR: ' sprintf('%.2f',PSNRs(ii)) '\n']);
 394         %     end
 395         %
 396
 397         fprintf('FINAL ESTIMATE, PSNR: %.3f dB, ISNR: %.3f dB\n', ...
 398              PSNR_FINAL_ESTIMATE, ISNR_FINAL_ESTIMATE);
 399
 400         figure, imshow(double(z(:,:,ceil(NumberOfFrames/2)))); % show the central frame
 401         title(sprintf('Noisy frame #%d',ceil(NumberOfFrames/2)));
 402
 403         figure, imshow(double(y_hat_wi(:,:,ceil(NumberOfFrames/2)))); % show the central frame
 404         title(sprintf('Denoised frame #%d',ceil(NumberOfFrames/2)));
 405
 406         text_psnr = sprintf('_PSNR%.2f', PSNR_FINAL_ESTIMATE);
 407     end
 408
 409     fprintf('The denoising took: %.1f sec (%.4f sec/frame). ', ...
 410         wiener_elapsed_time+estimate_elapsed_time, (wiener_elapsed_time+estimate_elapsed_time)/NumberOfFrames);
 411
 412
 413     text_vid = 'Denoised';
 414     FRATE = 30; % default value
 415     if isCharacterName,
 416         text_vid = Xnoisy_name(1:end-4);
 417         ainfo = aviinfo(Xnoisy_name);
 418         FRATE = ainfo.FramesPerSecond;
 419     end
 420
 421     avi_filename = sprintf('%s%s_%s_BM3D.avi', text_vid, text_psnr, bm3dProfile);
 422
 423     if exist(avi_filename, 'file') ~= 0,
 424         delete(avi_filename);
 425     end
 426     mov = avifile(avi_filename, 'Colormap', gray(256), 'compression', 'None', 'fps', FRATE);
 427     for ii = [1:NumberOfFrames],
 428         mov = addframe(mov, uint8(round(255*double(y_hat_wi(:,:,ii)))));
 429     end
 430     mov = close(mov);
 431     fprintf('The denoised video written to: %s.\n\n', avi_filename);
 432
 433 end
 434
 435 return;
 436
 437
 438
 439
 440
 441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 442 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 443 % Some auxiliary functions
 444 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 445 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 446
 447
 448
 449
 450 function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 451 %
 452 % Create forward and inverse transform matrices, which allow for perfect
 453 % reconstruction. The forward transform matrix is normalized so that the
 454 % l2-norm of each basis element is 1.
 455 %
 456 % [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 457 %
 458 %  INPUTS:
 459 %
 460 %   N               --> Size of the transform (for wavelets, must be 2^K)
 461 %
 462 %   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
 463 %                       listed by 'help wfilters' (bi-orthogonal wavelets)
 464 %                       'DCrand' -- an orthonormal transform with a DC and all
 465 %                       the other basis elements of random nature
 466 %
 467 %   dec_levels      --> If a wavelet transform is generated, this is the
 468 %                       desired decomposition level. Must be in the
 469 %                       range [0, log2(N)-1], where "0" implies
 470 %                       full decomposition.
 471 %
 472 %  OUTPUTS:
 473 %
 474 %   Tforward        --> (N x N) Forward transform matrix
 475 %
 476 %   Tinverse        --> (N x N) Inverse transform matrix
 477 %
 478
 479 if exist('dec_levels', 'var') ~= 1,
 480     dec_levels = 0;
 481 end
 482
 483 if N == 1,
 484     Tforward = 1;
 485 elseif strcmp(transform_type, 'hadamard') == 1,
 486     Tforward    = hadamard(N);
 487 elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
 488     Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 489        0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
 490        0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
 491       -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
 492        0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
 493                        0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
 494                        0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
 495                        0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
 496 elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 497     Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 498        0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
 499        0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
 500        0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
 501        0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
 502        0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
 503        0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
 504        0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
 505 elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
 506     Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
 507        0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
 508        0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
 509        0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
 510        0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
 511        0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
 512        0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
 513        0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
 514 elseif (N == 7) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 515     Tforward =[ 0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227   0.377964473009227;
 516        0.521120889169602   0.417906505941275   0.231920613924330                   0  -0.231920613924330  -0.417906505941275  -0.521120889169602;
 517        0.481588117120063   0.118942442321354  -0.333269317528993  -0.534522483824849  -0.333269317528993   0.118942442321354   0.481588117120063;
 518        0.417906505941275  -0.231920613924330  -0.521120889169602                   0   0.521120889169602   0.231920613924330  -0.417906505941275;
 519        0.333269317528993  -0.481588117120063  -0.118942442321354   0.534522483824849  -0.118942442321354  -0.481588117120063   0.333269317528993;
 520        0.231920613924330  -0.521120889169602   0.417906505941275                   0  -0.417906505941275   0.521120889169602  -0.231920613924330;
 521        0.118942442321354  -0.333269317528993   0.481588117120063  -0.534522483824849   0.481588117120063  -0.333269317528993   0.118942442321354];
 522 elseif strcmp(transform_type, 'dct') == 1,
 523     Tforward    = dct(eye(N));
 524 elseif strcmp(transform_type, 'dst') == 1,
 525     Tforward    = dst(eye(N));
 526 elseif strcmp(transform_type, 'DCrand') == 1,
 527     x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
 528     if (Q(1) < 0),
 529         Q = -Q;
 530     end;
 531     Tforward = Q';
 532 else %% a wavelet decomposition supported by 'wavedec'
 533     %%% Set periodic boundary conditions, to preserve bi-orthogonality
 534     dwtmode('per','nodisp');
 535
 536     Tforward = zeros(N,N);
 537     for i = 1:N
 538         Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
 539     end
 540 end
 541
 542 %%% Normalize the basis elements
 543 Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';
 544
 545 %%% Compute the inverse transform matrix
 546 Tinverse = inv(Tforward);
 547
 548 return;