BM3D/BM3D.m

   1 function [PSNR, y_est] = BM3D(y, z, sigma, profile, print_to_screen)
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %
   4 %  BM3D is an algorithm for attenuation of additive white Gaussian noise from
   5 %  grayscale images. This algorithm reproduces the results from the article:
   6 %
   7 %  [1] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Image Denoising
   8 %      by Sparse 3D Transform-Domain Collaborative Filtering,"
   9 %      IEEE Transactions on Image Processing, vol. 16, no. 8, August, 2007.
  10 %      preprint at http://www.cs.tut.fi/~foi/GCF-BM3D.
  11 %
  12 %
  13 %  FUNCTION INTERFACE:
  14 %
  15 %  [PSNR, y_est] = BM3D(y, z, sigma, profile, print_to_screen)
  16 %
  17 %  ! The function can work without any of the input arguments,
  18 %   in which case, the internal default ones are used !
  19 %
  20 %  BASIC USAGE EXAMPLES:
  21 %
  22 %     Case 1) Using the default parameters (i.e., image name, sigma, etc.)
  23 %
  24 %      [PSNR, y_est] = BM3D;
  25 %
  26 %     Case 2) Using an external noisy image:
  27 %
  28 %      % Read a grayscale image and scale its intensities in range [0,1]
  29 %      y = im2double(imread('Cameraman256.png'));
  30 %      % Generate the same seed used in the experimental results of [1]
  31 %      randn('seed', 0);
  32 %      % Standard deviation of the noise --- corresponding to intensity
  33 %      %  range [0,255], despite that the input was scaled in [0,1]
  34 %      sigma = 25;
  35 %      % Add the AWGN with zero mean and standard deviation 'sigma'
  36 %      z = y + (sigma/255)*randn(size(y));
  37 %      % Denoise 'z'. The denoised image is 'y_est', and 'NA = 1' because
  38 %      %  the true image was not provided
  39 %      [NA, y_est] = BM3D(1, z, sigma);
  40 %      % Compute the putput PSNR
  41 %      PSNR = 10*log10(1/mean((y(:)-y_est(:)).^2))
  42 %      % show the noisy image 'z' and the denoised 'y_est'
  43 %      figure; imshow(z);
  44 %      figure; imshow(y_est);
  45 %
  46 %     Case 3) If the original image y is provided as the first input
  47 %      argument, then some additional information is printed (PSNRs,
  48 %      figures, etc.). That is, "[NA, y_est] = BM3D(1, z, sigma);" in the
  49 %      above code should be replaced with:
  50 %
  51 %      [PSNR, y_est] = BM3D(y, z, sigma);
  52 %
  53 %
  54 %  INPUT ARGUMENTS (OPTIONAL):
  55 %
  56 %     1) y (matrix M x N): Noise-free image (needed for computing PSNR),
  57 %                           replace with the scalar 1 if not available.
  58 %     2) z (matrix M x N): Noisy image (intensities in range [0,1] or [0,255])
  59 %     3) sigma (double)  : Std. dev. of the noise (corresponding to intensities
  60 %                          in range [0,255] even if the range of z is [0,1])
  61 %     4) profile (char)  : 'np' --> Normal Profile
  62 %                          'lc' --> Fast Profile
  63 %     5) print_to_screen : 0 --> do not print output information (and do
  64 %                                not plot figures)
  65 %                          1 --> print information and plot figures
  66 %
  67 %  OUTPUTS:
  68 %     1) PSNR (double)          : Output PSNR (dB), only if the original
  69 %                                 image is available, otherwise PSNR = 0
  70 %     2) y_est (matrix M x N): Final estimate (in the range [0,1])
  71 %
  72 %
  73 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  74 %
  75 % Copyright (c) 2006-2011 Tampere University of Technology.
  76 % All rights reserved.
  77 % This work should only be used for nonprofit purposes.
  78 %
  79 % AUTHORS:
  80 %     Kostadin Dabov, email: dabov _at_ cs.tut.fi
  81 %
  82 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  83
  84
  85 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  86 %%%% In case, a noisy image z is not provided, then use the filename
  87 %%%%  below to read an original image (might contain path also). Later,
  88 %%%%  artificial AWGN noise is added and this noisy image is processed
  89 %%%%  by the BM3D.
  90 %%%%
  91 image_name = [
  92 %     'montage.png'
  93      'Cameraman256.png'
  94 %     'boat.png'
  95 %     'Lena512.png'
  96 %     'house.png'
  97 %     'barbara.png'
  98 %     'peppers256.png'
  99 %     'fingerprint.png'
 100 %     'couple.png'
 101 %     'hill.png'
 102 %     'man.png'
 103     ];
 104
 105 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 106 %%%%  Quality/complexity trade-off profile selection
 107 %%%%
 108 %%%%  'np' --> Normal Profile (balanced quality)
 109 %%%%  'lc' --> Low Complexity Profile (fast, lower quality)
 110 %%%%
 111 %%%%  'high' --> High Profile (high quality, not documented in [1])
 112 %%%%
 113 %%%%  'vn' --> This profile is automatically enabled for high noise
 114 %%%%           when sigma > 40
 115 %%%%
 116 %%%%  'vn_old' --> This is the old 'vn' profile that was used in [1].
 117 %%%%           It gives inferior results than 'vn' in most cases.
 118 %%%%
 119 if (exist('profile') ~= 1)
 120     profile         = 'np'; %% default profile
 121 end
 122
 123 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 124 %%%%  Specify the std. dev. of the corrupting noise
 125 %%%%
 126 if (exist('sigma') ~= 1),
 127     sigma               = 25; %% default standard deviation of the AWGN
 128 end
 129
 130 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 131 %%%% Following are the parameters for the Normal Profile.
 132 %%%%
 133
 134 %%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
 135 transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
 136 transform_2D_Wiener_name = 'dct';     %% transform used for the Wiener filt. of size N1_wiener x N1_wiener
 137 transform_3rd_dim_name   = 'haar';    %% transform used in the 3-rd dim, the same for HT and Wiener filt.
 138
 139 %%%% Hard-thresholding (HT) parameters:
 140 N1                  = 8;   %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
 141 Nstep               = 3;   %% sliding step to process every next reference block
 142 N2                  = 16;  %% maximum number of similar blocks (maximum size of the 3rd dimension of a 3D array)
 143 Ns                  = 39;  %% length of the side of the search neighborhood for full-search block-matching (BM), must be odd
 144 tau_match           = 3000;%% threshold for the block-distance (d-distance)
 145 lambda_thr2D        = 0;   %% threshold parameter for the coarse initial denoising used in the d-distance measure
 146 lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D transform domain
 147 beta                = 2.0; %% parameter of the 2D Kaiser window used in the reconstruction
 148
 149 %%%% Wiener filtering parameters:
 150 N1_wiener           = 8;
 151 Nstep_wiener        = 3;
 152 N2_wiener           = 32;
 153 Ns_wiener           = 39;
 154 tau_match_wiener    = 400;
 155 beta_wiener         = 2.0;
 156
 157 %%%% Block-matching parameters:
 158 stepFS              = 1;  %% step that forces to switch to full-search BM, "1" implies always full-search
 159 smallLN             = 'not used in np'; %% if stepFS > 1, then this specifies the size of the small local search neighb.
 160 stepFSW             = 1;
 161 smallLNW            = 'not used in np';
 162 thrToIncStep        = 8;  % if the number of non-zero coefficients after HT is less than thrToIncStep,
 163                           % than the sliding step to the next reference block is incresed to (nm1-1)
 164
 165 if strcmp(profile, 'lc') == 1,
 166
 167     Nstep               = 6;
 168     Ns                  = 25;
 169     Nstep_wiener        = 5;
 170     N2_wiener           = 16;
 171     Ns_wiener           = 25;
 172
 173     thrToIncStep        = 3;
 174     smallLN             = 3;
 175     stepFS              = 6*Nstep;
 176     smallLNW            = 2;
 177     stepFSW             = 5*Nstep_wiener;
 178
 179 end
 180
 181 % Profile 'vn' was proposed in
 182 %  Y. Hou, C. Zhao, D. Yang, and Y. Cheng, 'Comment on "Image Denoising by Sparse 3D Transform-Domain
 183 %  Collaborative Filtering"', accepted for publication, IEEE Trans. on Image Processing, July, 2010.
 184 % as a better alternative to that initially proposed in [1] (which is currently in profile 'vn_old')
 185 if (strcmp(profile, 'vn') == 1) | (sigma > 40),
 186
 187     N2                  = 32;
 188     Nstep               = 4;
 189
 190     N1_wiener           = 11;
 191     Nstep_wiener        = 6;
 192
 193     lambda_thr3D        = 2.8;
 194     thrToIncStep        = 3;
 195     tau_match_wiener    = 3500;
 196     tau_match           = 25000;
 197
 198     Ns_wiener           = 39;
 199
 200 end
 201
 202 % The 'vn_old' profile corresponds to the original parameters for strong noise proposed in [1].
 203 if (strcmp(profile, 'vn_old') == 1) & (sigma > 40),
 204
 205     transform_2D_HT_name = 'dct';
 206
 207     N1                  = 12;
 208     Nstep               = 4;
 209
 210     N1_wiener           = 11;
 211     Nstep_wiener        = 6;
 212
 213     lambda_thr3D        = 2.8;
 214     lambda_thr2D        = 2.0;
 215     thrToIncStep        = 3;
 216     tau_match_wiener    = 3500;
 217     tau_match           = 5000;
 218
 219     Ns_wiener           = 39;
 220
 221 end
 222
 223 decLevel = 0;        %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
 224 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
 225
 226 if strcmp(profile, 'high') == 1, %% this profile is not documented in [1]
 227
 228     decLevel     = 1;
 229     Nstep        = 2;
 230     Nstep_wiener = 2;
 231     lambda_thr3D = 2.5;
 232     vMask = ones(N1,1); vMask((end/4+1):end/2)= 1.01; vMask((end/2+1):end) = 1.07; %% this allows to have different threhsolds for the finest and next-to-the-finest subbands
 233     thr_mask = vMask * vMask';
 234     beta         = 2.5;
 235     beta_wiener  = 1.5;
 236
 237 end
 238
 239 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 240 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 241 %%%% Note: touch below this point only if you know what you are doing!
 242 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 243
 244 %%% Check whether to dump information to the screen or remain silent
 245 dump_output_information = 1;
 246 if (exist('print_to_screen') == 1) & (print_to_screen == 0),
 247     dump_output_information = 0;
 248 end
 249
 250 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 251 %%%% Create transform matrices, etc.
 252 %%%%
 253 [Tfor, Tinv]   = getTransfMatrix(N1, transform_2D_HT_name, decLevel);     %% get (normalized) forward and inverse transform matrices
 254 [TforW, TinvW] = getTransfMatrix(N1_wiener, transform_2D_Wiener_name, 0); %% get (normalized) forward and inverse transform matrices
 255
 256 if (strcmp(transform_3rd_dim_name, 'haar') == 1) | (strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
 257     %%% If Haar is used in the 3-rd dimension, then a fast internal transform is used, thus no need to generate transform
 258     %%% matrices.
 259     hadper_trans_single_den         = {};
 260     inverse_hadper_trans_single_den = {};
 261 else
 262     %%% Create transform matrices. The transforms are later applied by
 263     %%% matrix-vector multiplication for the 1D case.
 264     for hpow = 0:ceil(log2(max(N2,N2_wiener))),
 265         h = 2^hpow;
 266         [Tfor3rd, Tinv3rd]   = getTransfMatrix(h, transform_3rd_dim_name, 0);
 267         hadper_trans_single_den{h}         = single(Tfor3rd);
 268         inverse_hadper_trans_single_den{h} = single(Tinv3rd');
 269     end
 270 end
 271
 272 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 273 %%%% 2D Kaiser windows used in the aggregation of block-wise estimates
 274 %%%%
 275 if beta_wiener==2 & beta==2 & N1_wiener==8 & N1==8 % hardcode the window function so that the signal processing toolbox is not needed by default
 276     Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
 277         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 278         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 279         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 280         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 281         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 282         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 283         0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924];
 284     Wwin2D_wiener = Wwin2D;
 285 else
 286     Wwin2D           = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
 287     Wwin2D_wiener    = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the aggregation of the Wiener filt. part
 288 end
 289 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 290 %%%% If needed, read images, generate noise, or scale the images to the
 291 %%%% [0,1] interval
 292 %%%%
 293 if (exist('y') ~= 1) | (exist('z') ~= 1)
 294     y        = im2double(imread(image_name));  %% read a noise-free image and put in intensity range [0,1]
 295     randn('seed', 0);                          %% generate seed
 296     z        = y + (sigma/255)*randn(size(y)); %% create a noisy image
 297 else  % external images
 298
 299     image_name = 'External image';
 300
 301     % convert z to double precision if needed
 302     z = double(z);
 303
 304     % convert y to double precision if needed
 305     y = double(y);
 306
 307     % if z's range is [0, 255], then convert to [0, 1]
 308     if (max(z(:)) > 10), % a naive check for intensity range
 309         z = z / 255;
 310     end
 311
 312     % if y's range is [0, 255], then convert to [0, 1]
 313     if (max(y(:)) > 10), % a naive check for intensity range
 314         y = y / 255;
 315     end
 316 end
 317
 318
 319
 320 if (size(z,3) ~= 1) | (size(y,3) ~= 1),
 321     error('BM3D accepts only grayscale 2D images.');
 322 end
 323
 324
 325 % Check if the true image y is a valid one; if not, then we cannot compute PSNR, etc.
 326 y_is_invalid_image = (length(size(z)) ~= length(size(y))) | (size(z,1) ~= size(y,1)) | (size(z,2) ~= size(y,2));
 327 if (y_is_invalid_image),
 328     dump_output_information = 0;
 329 end
 330
 331 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 332 %%% Print image information to the screen
 333 %%%%
 334 if dump_output_information == 1,
 335     fprintf('Image: %s (%dx%d), sigma: %.1f\n', image_name, size(z,1), size(z,2), sigma);
 336 end
 337
 338 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 339 %%%% Step 1. Produce the basic estimate by HT filtering
 340 %%%%
 341 tic;
 342 y_hat = bm3d_thr(z, hadper_trans_single_den, Nstep, N1, N2, lambda_thr2D,...
 343         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), Wwin2D, smallLN, stepFS );
 344 estimate_elapsed_time = toc;
 345
 346 if dump_output_information == 1,
 347     PSNR_INITIAL_ESTIMATE = 10*log10(1/mean((y(:)-double(y_hat(:))).^2));
 348     fprintf('BASIC ESTIMATE, PSNR: %.2f dB\n', PSNR_INITIAL_ESTIMATE);
 349 end
 350
 351
 352 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 353 %%%% Step 2. Produce the final estimate by Wiener filtering (using the
 354 %%%%  hard-thresholding initial estimate)
 355 %%%%
 356 tic;
 357 y_est = bm3d_wiener(z, y_hat, hadper_trans_single_den, Nstep_wiener, N1_wiener, N2_wiener, ...
 358     '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, Wwin2D_wiener, smallLNW, stepFSW, single(ones(N1_wiener)) );
 359 wiener_elapsed_time = toc;
 360
 361 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 362 %%%% Calculate the final estimate's PSNR, print it, and show the
 363 %%%% denoised image next to the noisy one
 364 %%%%
 365 y_est = double(y_est);
 366
 367 PSNR = 0; %% Remains 0 if the true image y is not available
 368 if (~y_is_invalid_image), % checks if y is a valid image
 369     PSNR = 10*log10(1/mean((y(:)-y_est(:)).^2)); % y is valid
 370 end
 371
 372 if dump_output_information == 1,
 373     fprintf('FINAL ESTIMATE (total time: %.1f sec), PSNR: %.2f dB\n', ...
 374         wiener_elapsed_time + estimate_elapsed_time, PSNR);
 375
 376     figure, imshow(z); title(sprintf('Noisy %s, PSNR: %.3f dB (sigma: %d)', ...
 377         image_name(1:end-4), 10*log10(1/mean((y(:)-z(:)).^2)), sigma));
 378
 379     figure, imshow(y_est); title(sprintf('Denoised %s, PSNR: %.3f dB', ...
 380         image_name(1:end-4), PSNR));
 381
 382 end
 383
 384 return;
 385
 386
 387
 388
 389
 390 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 391 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 392 % Some auxiliary functions
 393 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 394 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 395
 396
 397
 398
 399 function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 400 %
 401 % Create forward and inverse transform matrices, which allow for perfect
 402 % reconstruction. The forward transform matrix is normalized so that the
 403 % l2-norm of each basis element is 1.
 404 %
 405 % [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 406 %
 407 %  INPUTS:
 408 %
 409 %   N               --> Size of the transform (for wavelets, must be 2^K)
 410 %
 411 %   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
 412 %                       listed by 'help wfilters' (bi-orthogonal wavelets)
 413 %                       'DCrand' -- an orthonormal transform with a DC and all
 414 %                       the other basis elements of random nature
 415 %
 416 %   dec_levels      --> If a wavelet transform is generated, this is the
 417 %                       desired decomposition level. Must be in the
 418 %                       range [0, log2(N)-1], where "0" implies
 419 %                       full decomposition.
 420 %
 421 %  OUTPUTS:
 422 %
 423 %   Tforward        --> (N x N) Forward transform matrix
 424 %
 425 %   Tinverse        --> (N x N) Inverse transform matrix
 426 %
 427
 428 if exist('dec_levels') ~= 1,
 429     dec_levels = 0;
 430 end
 431
 432 if N == 1,
 433     Tforward = 1;
 434 elseif strcmp(transform_type, 'hadamard') == 1,
 435     Tforward    = hadamard(N);
 436 elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
 437     Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 438        0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
 439        0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
 440       -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
 441        0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
 442                        0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
 443                        0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
 444                        0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
 445 elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 446     Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 447        0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
 448        0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
 449        0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
 450        0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
 451        0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
 452        0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
 453        0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
 454 elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
 455     Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
 456        0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
 457        0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
 458        0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
 459        0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
 460        0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
 461        0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
 462        0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
 463 elseif strcmp(transform_type, 'dct') == 1,
 464     Tforward    = dct(eye(N));
 465 elseif strcmp(transform_type, 'dst') == 1,
 466     Tforward    = dst(eye(N));
 467 elseif strcmp(transform_type, 'DCrand') == 1,
 468     x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
 469     if (Q(1) < 0),
 470         Q = -Q;
 471     end;
 472     Tforward = Q';
 473 else %% a wavelet decomposition supported by 'wavedec'
 474     %%% Set periodic boundary conditions, to preserve bi-orthogonality
 475     dwtmode('per','nodisp');
 476
 477     Tforward = zeros(N,N);
 478     for i = 1:N
 479         Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
 480     end
 481 end
 482
 483 %%% Normalize the basis elements
 484 Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';
 485
 486 %%% Compute the inverse transform matrix
 487 Tinverse = inv(Tforward);
 488
 489 return;
 490