BM3D/CBM3D.m

   1 function [PSNR, yRGB_est] = CBM3D(yRGB, zRGB, sigma, profile, print_to_screen, colorspace)
   2 %
   3 %  CBM3D is algorithm for attenuation of additive white Gaussian noise from
   4 %  color RGB images. This algorithm reproduces the results from the article:
   5 %
   6 %  [1] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Color image
   7 %   denoising via sparse 3D collaborative filtering with grouping constraint in
   8 %   luminance-chrominance space," submitted to IEEE Int. Conf. Image Process.,
   9 %   January 2007, in review, preprint at http://www.cs.tut.fi/~foi/GCF-BM3D.
  10 %
  11 %  FUNCTION INTERFACE:
  12 %
  13 %  [PSNR, yRGB_est] = CBM3D(yRGB, zRGB, sigma, profile, print_to_screen, colorspace)
  14 %
  15 %  ! The function can work without any of the input arguments,
  16 %   in which case, the internal default ones are used !
  17 %
  18 %  BASIC USAGE EXAMPLES:
  19 %
  20 %     Case 1) Using the default parameters (i.e., image name, sigma, etc.)
  21 %
  22 %      [PSNR, yRGB_est] = CBM3D;
  23 %
  24 %     Case 2) Using an external noisy image:
  25 %
  26 %      % Read an RGB image and scale its intensities in range [0,1]
  27 %      yRGB = im2double(imread('image_House256rgb.png'));
  28 %      % Generate the same seed used in the experimental results of [1]
  29 %      randn('seed', 0);
  30 %      % Standard deviation of the noise --- corresponding to intensity
  31 %      %  range [0,255], despite that the input was scaled in [0,1]
  32 %      sigma = 25;
  33 %      % Add the AWGN with zero mean and standard deviation 'sigma'
  34 %      zRGB = yRGB + (sigma/255)*randn(size(yRGB));
  35 %      % Denoise 'zRGB'. The denoised image is 'yRGB_est', and 'NA = 1'
  36 %      %  because the true image was not provided
  37 %      [NA, yRGB_est] = CBM3D(1, zRGB, sigma);
  38 %      % Compute the putput PSNR
  39 %      PSNR = 10*log10(1/mean((yRGB(:)-yRGB_est(:)).^2))
  40 %      % show the noisy image 'zRGB' and the denoised 'yRGB_est'
  41 %      figure; imshow(min(max(zRGB,0),1));
  42 %      figure; imshow(min(max(yRGB_est,0),1));
  43 %
  44 %     Case 3) If the original image yRGB is provided as the first input
  45 %      argument, then some additional information is printed (PSNRs,
  46 %      figures, etc.). That is, "[NA, yRGB_est] = BM3D(1, zRGB, sigma);" in the
  47 %      above code should be replaced with:
  48 %
  49 %      [PSNR, yRGB_est] = CBM3D(yRGB, zRGB, sigma);
  50 %
  51 %
  52 %  INPUT ARGUMENTS (OPTIONAL):
  53 %     1) yRGB (M x N x 3): Noise-free RGB image (needed for computing PSNR),
  54 %                           replace with the scalar 1 if not available.
  55 %     2) zRGB (M x N x 3): Noisy RGBimage (intensities in range [0,1] or [0,255])
  56 %     3) sigma (double)  : Std. dev. of the noise (corresponding to intensities
  57 %                            in range [0,255] even if the range of zRGB is [0,1])
  58 %     4) profile (char)  : 'np' --> Normal Profile
  59 %                          'lc' --> Fast Profile
  60 %     5) print_to_screen : 0 --> do not print output information (and do
  61 %                                not plot figures)
  62 %                          1 --> print information and plot figures
  63 %     6) colorspace (char): 'opp'   --> use opponent colorspace
  64 %                          'yCbCr' --> use yCbCr colorspace
  65 %
  66 %  OUTPUTS:
  67 %     1) PSNR (double)          : Output PSNR (dB), only if the original
  68 %                                 image is available, otherwise PSNR = 0
  69 %     2) yRGB_est (M x N x 3): Final RGB estimate (in the range [0,1])
  70 %
  71 %
  72 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  73 %
  74 % Copyright (c) 2007-2011 Tampere University of Technology.
  75 % All rights reserved.
  76 % This work should only be used for nonprofit purposes.
  77 %
  78 % AUTHORS:
  79 %     Kostadin Dabov, email: dabov _at_ cs.tut.fi
  80 %
  81 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  82
  83 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  84 %%%% In case, there is no input image (zRGB or yRGB), then use the filename
  85 %%%%  below to read an original image (might contain path also). Later,
  86 %%%%  artificial AWGN noise is added and this noisy image is processed
  87 %%%%  by the CBM3D.
  88 %%%%
  89 image_name = [
  90 %    'kodim12.png'
  91    'image_Lena512rgb.png'
  92 %     'image_House256rgb.png'
  93 %    'image_Peppers512rgb.png'
  94 %    'image_Baboon512rgb.png'
  95 %    'image_F16_512rgb.png'
  96    ];
  97
  98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  99 %%%%  Quality/complexity trade-off
 100 %%%%
 101 %%%%  'np' --> Normal Profile (balanced quality)
 102 %%%%  'lc' --> Low Complexity Profile (fast, lower quality)
 103 %%%%
 104 %%%%  'high' --> High Profile (high quality, not documented in [1])
 105 %%%%
 106 %%%%  'vn' --> This profile is automatically enabled for high noise
 107 %%%%           when sigma > 40
 108 %%%%
 109 %%%%  'vn_old' --> This is the old 'vn' profile that was used in [1].
 110 %%%%           It gives inferior results than 'vn' in most cases.
 111 %%%%
 112 if (exist('profile') ~= 1)
 113     profile         = 'np'; %% default profile
 114 end
 115
 116
 117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 118 %%%%  Specify the std. dev. of the corrupting noise
 119 %%%%
 120 if (exist('sigma') ~= 1),
 121    sigma                = 50; %% default standard deviation of the AWGN
 122 end
 123
 124 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 125 %%%%  Colorspace in which we perform denoising. BM is applied to the first
 126 %%%%  component and the matching information is re-used for the other two.
 127 %%%%
 128 if (exist('colorspace') ~= 1),
 129     colorspace              = 'opp'; %%% (valid colorspaces are: 'yCbCr' and 'opp')
 130 end
 131
 132 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 133 %%%% Following are the parameters for the Normal Profile.
 134 %%%%
 135
 136 %%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
 137 transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
 138 transform_2D_Wiener_name = 'dct';     %% transform used for the Wiener filt. of size N1_wiener x N1_wiener
 139 transform_3rd_dim_name   = 'haar';    %% transform used in the 3-rd dim, the same for HT and Wiener filt.
 140
 141 %%%% Hard-thresholding (HT) parameters:
 142 N1                  = 8;   %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
 143 Nstep               = 3;   %% sliding step to process every next reference block
 144 N2                  = 16;  %% maximum number of similar blocks (maximum size of the 3rd dimension of a 3D array)
 145 Ns                  = 39;  %% length of the side of the search neighborhood for full-search block-matching (BM), must be odd
 146 tau_match           = 3000;%% threshold for the block-distance (d-distance)
 147 lambda_thr2D        = 0;   %% threshold parameter for the coarse initial denoising used in the d-distance measure
 148 lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D transform domain
 149 beta                = 2.0; %% parameter of the 2D Kaiser window used in the reconstruction
 150
 151 %%%% Wiener filtering parameters:
 152 N1_wiener           = 8;
 153 Nstep_wiener        = 3;
 154 N2_wiener           = 32;
 155 Ns_wiener           = 39;
 156 tau_match_wiener    = 400;
 157 beta_wiener         = 2.0;
 158
 159 %%%% Block-matching parameters:
 160 stepFS              = 1;  %% step that forces to switch to full-search BM, "1" implies always full-search
 161 smallLN             = 'not used in np'; %% if stepFS > 1, then this specifies the size of the small local search neighb.
 162 stepFSW             = 1;
 163 smallLNW            = 'not used in np';
 164 thrToIncStep        = 8;  %% used in the HT filtering to increase the sliding step in uniform regions
 165
 166 if strcmp(profile, 'lc') == 1,
 167
 168     Nstep               = 6;
 169     Ns                  = 25;
 170     Nstep_wiener        = 5;
 171     N2_wiener           = 16;
 172     Ns_wiener           = 25;
 173
 174     thrToIncStep        = 3;
 175     smallLN             = 3;
 176     stepFS              = 6*Nstep;
 177     smallLNW            = 2;
 178     stepFSW             = 5*Nstep_wiener;
 179
 180 end
 181
 182 % Profile 'vn' was proposed in
 183 %  Y. Hou, C. Zhao, D. Yang, and Y. Cheng, 'Comment on "Image Denoising by Sparse 3D Transform-Domain
 184 %  Collaborative Filtering"', accepted for publication, IEEE Trans. on Image Processing, July, 2010.
 185 % as a better alternative to that initially proposed in [1] (which is currently in profile 'vn_old')
 186 if (strcmp(profile, 'vn') == 1) | (sigma > 40),
 187
 188     N2                  = 32;
 189     Nstep               = 4;
 190
 191     N1_wiener           = 11;
 192     Nstep_wiener        = 6;
 193
 194     lambda_thr3D        = 2.8;
 195     thrToIncStep        = 3;
 196     tau_match_wiener    = 3500;
 197     tau_match           = 25000;
 198
 199     Ns_wiener           = 39;
 200
 201 end
 202
 203 % The 'vn_old' profile corresponds to the original parameters for strong noise proposed in [1].
 204 if (strcmp(profile, 'vn_old') == 1) & (sigma > 40),
 205
 206     transform_2D_HT_name = 'dct';
 207
 208     N1                  = 12;
 209     Nstep               = 4;
 210
 211     N1_wiener           = 11;
 212     Nstep_wiener        = 6;
 213
 214     lambda_thr3D        = 2.8;
 215     lambda_thr2D        = 2.0;
 216     thrToIncStep        = 3;
 217     tau_match_wiener    = 3500;
 218     tau_match           = 5000;
 219
 220     Ns_wiener           = 39;
 221
 222 end
 223
 224
 225 decLevel = 0;        %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
 226 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
 227
 228 if strcmp(profile, 'high') == 1,
 229
 230     decLevel     = 1;
 231     Nstep        = 2;
 232     Nstep_wiener = 2;
 233     lambda_thr3D = 2.5;
 234     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
 235     thr_mask = vMask * vMask';
 236     beta         = 2.5;
 237     beta_wiener  = 1.5;
 238
 239 end
 240
 241 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 242 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 243 %%%% Note: touch below this point only if you know what you are doing!
 244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 245 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 246
 247 %%% Check whether to dump information to the screen or reamin silent
 248 dump_output_information = 1;
 249 if (exist('print_to_screen') == 1) & (print_to_screen == 0),
 250     dump_output_information = 0;
 251 end
 252
 253 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 254 %%%% Create transform matrices, etc.
 255 %%%%
 256 [Tfor, Tinv]   = getTransfMatrix(N1, transform_2D_HT_name, decLevel);  %% get (normalized) forward and inverse transform matrices
 257 [TforW, TinvW] = getTransfMatrix(N1_wiener, transform_2D_Wiener_name); %% get (normalized) forward and inverse transform matrices
 258
 259 if (strcmp(transform_3rd_dim_name, 'haar') == 1) | (strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
 260     %%% If Haar is used in the 3-rd dimension, then a fast internal transform is used, thus no need to generate transform
 261     %%% matrices.
 262     hadper_trans_single_den         = {};
 263     inverse_hadper_trans_single_den = {};
 264 else
 265     %%% Create transform matrices. The transforms are later applied by
 266     %%% matrix-vector multiplication for the 1D case.
 267     for hpow = 0:ceil(log2(max(N2,N2_wiener))),
 268         h = 2^hpow;
 269         [Tfor3rd, Tinv3rd]   = getTransfMatrix(h, transform_3rd_dim_name, 0);
 270         hadper_trans_single_den{h}         = single(Tfor3rd);
 271         inverse_hadper_trans_single_den{h} = single(Tinv3rd');
 272     end
 273 end
 274
 275 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 276 %%%% 2D Kaiser windows used in the aggregation of block-wise estimates
 277 %%%%
 278 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
 279     Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
 280         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 281         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 282         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 283         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 284         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 285         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 286         0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924];
 287     Wwin2D_wiener = Wwin2D;
 288 else
 289     Wwin2D           = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
 290     Wwin2D_wiener    = kaiser(N1_wiener, beta_wiener) * kaiser(N1_wiener, beta_wiener)'; % Kaiser window used in the aggregation of the Wiener filt. part
 291 end
 292
 293 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 294 %%%% If needed, read images, generate noise, or scale the images to the
 295 %%%% [0,1] interval
 296 %%%%
 297 if (exist('yRGB') ~= 1) | (exist('zRGB') ~= 1)
 298     yRGB        = im2double(imread(image_name));  %% read a noise-free image
 299     randn('seed', 0);                          %% generate seed
 300     zRGB        = yRGB + (sigma/255)*randn(size(yRGB)); %% create a noisy image
 301 else % external images
 302     image_name = 'External image';
 303
 304     % convert zRGB to double precision
 305     zRGB = double(zRGB);
 306
 307     % convert yRGB to double precision
 308     yRGB = double(yRGB);
 309
 310     % if zRGB's range is [0, 255], then convert to [0, 1]
 311     if (max(zRGB(:)) > 10), % a naive check for intensity range
 312         zRGB = zRGB / 255;
 313     end
 314
 315     % if yRGB's range is [0, 255], then convert to [0, 1]
 316     if (max(yRGB(:)) > 10), % a naive check for intensity range
 317         yRGB = yRGB / 255;
 318     end
 319 end
 320
 321
 322 if (size(zRGB,3) ~= 3) | (size(zRGB,4) ~= 1),
 323     error('CBM3D accepts only input RGB images (i.e. matrices of size M x N x 3).');
 324 end
 325
 326 % Check if the true image yRGB is a valid one; if not, then we cannot compute PSNR, etc.
 327 yRGB_is_invalid_image = (length(size(zRGB)) ~= length(size(yRGB))) | (size(zRGB,1) ~= size(yRGB,1)) | (size(zRGB,2) ~= size(yRGB,2)) | (size(zRGB,3) ~= size(yRGB,3));
 328 if (yRGB_is_invalid_image),
 329     dump_output_information = 0;
 330 end
 331
 332
 333 [Xv, Xh, numSlices] = size(zRGB);              %%% obtain image sizes
 334
 335 if numSlices ~= 3
 336     fprintf('Error, an RGB color image is required!\n');
 337     return;
 338 end
 339
 340 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 341 %%%% Change colorspace, compute the l2-norms of the new color channels
 342 %%%%
 343 [zColSpace l2normLumChrom] = function_rgb2LumChrom(zRGB, colorspace);
 344
 345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 346 %%%% Print image information to the screen
 347 %%%%
 348 if dump_output_information == 1,
 349     fprintf(sprintf('Image: %s (%dx%dx%d), sigma: %.1f\n', image_name, Xv, Xh, numSlices, sigma));
 350 end
 351
 352
 353 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 354 %%%% Step 1. Basic estimate by collaborative hard-thresholding and using
 355 %%%% the grouping constraint on the chrominances.
 356 %%%%
 357 tic;
 358 y_hat = bm3d_thr_color(zColSpace, hadper_trans_single_den, Nstep, N1, N2, lambda_thr2D,...
 359     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', 'unused arg', l2normLumChrom, Wwin2D, smallLN, stepFS );
 360 estimate_elapsed_time = toc;
 361
 362
 363 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 364 %%%% Step 2. Final estimate by collaborative Wiener filtering and using
 365 %%%% the grouping constraint on the chrominances.
 366 %%%%
 367 tic;
 368 yRGB_est = bm3d_wiener_color(zColSpace, y_hat, hadper_trans_single_den, Nstep_wiener, N1_wiener, N2_wiener, ...
 369     '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', 'unused arg', l2normLumChrom, Wwin2D_wiener, smallLNW, stepFSW );
 370 wiener_elapsed_time = toc;
 371
 372 yRGB_est = double(yRGB_est);
 373
 374 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 375 %%%% Convert back to RGB colorspace
 376 %%%%
 377 yRGB_est = function_LumChrom2rgb(yRGB_est, colorspace);
 378
 379
 380 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 381 %%%% Calculate final estimate's PSNR and ISNR, print them, and show the
 382 %%%% denoised image
 383 %%%%
 384 PSNR = 0; %% Remains 0 if the true image yRGB is not available
 385 if (~yRGB_is_invalid_image), % then we assume yRGB is a valid image
 386     PSNR = 10*log10(1/mean((yRGB(:)-yRGB_est(:)).^2));
 387 end
 388
 389 if dump_output_information == 1,
 390     fprintf(sprintf('FINAL ESTIMATE (total time: %.1f sec), PSNR: %.2f dB\n', ...
 391         wiener_elapsed_time + estimate_elapsed_time, PSNR));
 392
 393     figure, imshow(min(max(zRGB,0),1)); title(sprintf('Noisy %s, PSNR: %.3f dB (sigma: %d)', ...
 394         image_name(1:end-4), 10*log10(1/mean((yRGB(:)-zRGB(:)).^2)), sigma));
 395
 396     figure, imshow(min(max(yRGB_est,0),1)); title(sprintf('Denoised %s, PSNR: %.3f dB', ...
 397         image_name(1:end-4), PSNR));
 398 end
 399
 400 return;
 401
 402
 403
 404
 405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 406 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 407 % Some auxiliary functions
 408 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 409 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 410
 411
 412
 413
 414 function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 415 %
 416 % Create forward and inverse transform matrices, which allow for perfect
 417 % reconstruction. The forward transform matrix is normalized so that the
 418 % l2-norm of each basis element is 1.
 419 %
 420 % [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 421 %
 422 %  INPUTS:
 423 %
 424 %   N               --> Size of the transform (for wavelets, must be 2^K)
 425 %
 426 %   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
 427 %                       listed by 'help wfilters' (bi-orthogonal wavelets)
 428 %                       'DCrand' -- an orthonormal transform with a DC and all
 429 %                       the other basis elements of random nature
 430 %
 431 %   dec_levels      --> If a wavelet transform is generated, this is the
 432 %                       desired decomposition level. Must be in the
 433 %                       range [0, log2(N)-1], where "0" implies
 434 %                       full decomposition.
 435 %
 436 %  OUTPUTS:
 437 %
 438 %   Tforward        --> (N x N) Forward transform matrix
 439 %
 440 %   Tinverse        --> (N x N) Inverse transform matrix
 441 %
 442
 443 if exist('dec_levels') ~= 1,
 444     dec_levels = 0;
 445 end
 446
 447 if N == 1,
 448     Tforward = 1;
 449 elseif strcmp(transform_type, 'hadamard') == 1,
 450     Tforward    = hadamard(N);
 451 elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
 452     Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 453        0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
 454        0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
 455       -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
 456        0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
 457                        0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
 458                        0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
 459                        0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
 460 elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 461     Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 462        0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
 463        0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
 464        0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
 465        0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
 466        0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
 467        0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
 468        0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
 469 elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
 470     Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
 471        0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
 472        0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
 473        0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
 474        0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
 475        0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
 476        0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
 477        0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
 478 elseif strcmp(transform_type, 'dct') == 1,
 479     Tforward    = dct(eye(N));
 480 elseif strcmp(transform_type, 'dst') == 1,
 481     Tforward    = dst(eye(N));
 482 elseif strcmp(transform_type, 'DCrand') == 1,
 483     x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
 484     if (Q(1) < 0),
 485         Q = -Q;
 486     end;
 487     Tforward = Q';
 488 else %% a wavelet decomposition supported by 'wavedec'
 489     %%% Set periodic boundary conditions, to preserve bi-orthogonality
 490     dwtmode('per','nodisp');
 491
 492     Tforward = zeros(N,N);
 493     for i = 1:N
 494         Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
 495     end
 496 end
 497
 498 %%% Normalize the basis elements
 499 Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';
 500
 501 %%% Compute the inverse transform matrix
 502 Tinverse = inv(Tforward);
 503
 504 return;
 505
 506 function [y, A, l2normLumChrom]=function_rgb2LumChrom(xRGB, colormode)
 507 % Forward color-space transformation   ( inverse transformation is function_LumChrom2rgb.m )
 508 %
 509 % Alessandro Foi - Tampere University of Technology - 2005 - 2006   Public release v1.03 (March 2006)
 510 % -----------------------------------------------------------------------------------------------------------------------------------------------
 511 %
 512 % SYNTAX:
 513 %
 514 %   [y A l2normLumChrom] = function_rgb2LumChrom(xRGB, colormode);
 515 %
 516 % INPUTS:
 517 %   xRGB  is RGB image with range [0 1]^3
 518 %
 519 %   colormode = 'opp', 'yCbCr', 'pca', or a custom 3x3 matrix
 520 %
 521 %       'opp'     Opponent color space ('opp' is equirange version)
 522 %       'yCbCr'   The standard yCbCr (e.g. for JPEG images)
 523 %       'pca'     Principal components   (note that this transformation is renormalized to be equirange)
 524 %
 525 % OUTPUTS:
 526 %   y  is color-transformed image (with range typically included in or equal to [0 1]^3, depending on the transformation matrix)
 527 %
 528 %   l2normLumChrom (optional) l2-norm of the transformation (useful for noise std calculation)
 529 %   A  transformation matrix  (used necessarily if colormode='pca')
 530 %
 531 %   NOTES:  -  If only two outputs are used, then the second output is l2normLumChrom, unless colormode='pca';
 532 %           -  'opp' is used by default if no colormode is specified.
 533 %
 534 %
 535 % USAGE EXAMPLE FOR PCA TRANSFORMATION:
 536 %  %%%%  -- forward color transformation --
 537 %    if colormode=='pca'
 538 %       [zLumChrom colormode] = function_rgb2LumChrom(zRGB,colormode); % 'colormode' is assigned a 3x3 transform matrix
 539 %    else
 540 %       zLumChrom = function_rgb2LumChrom(zRGB,colormode);
 541 %    end
 542 %
 543 %  %%%% [ ... ]  Some processing  [ ... ]
 544 %
 545 %  %%%%  -- inverse color transformation --
 546 %    zRGB=function_LumChrom2rgb(zLumChrom,colormode);
 547 %
 548
 549 if nargin==1
 550     colormode='opp';
 551 end
 552 change_output=0;
 553 if size(colormode)==[3 3]
 554     A=colormode;
 555     l2normLumChrom=sqrt(sum(A.^2,2));
 556 else
 557     if strcmp(colormode,'opp')
 558         A=[1/3 1/3 1/3; 0.5  0  -0.5; 0.25  -0.5  0.25];
 559     end
 560     if strcmp(colormode,'yCbCr')
 561         A=[0.299   0.587   0.114;   -0.16873660714285  -0.33126339285715   0.5;   0.5  -0.4186875  -0.0813125];
 562     end
 563     if strcmp(colormode,'pca')
 564         A=princomp(reshape(xRGB,[size(xRGB,1)*size(xRGB,2) 3]))';
 565         A=A./repmat(sum(A.*(A>0),2)-sum(A.*(A<0),2),[1 3]);  %% ranges are normalized to unitary length;
 566     else
 567         if nargout==2
 568             change_output=1;
 569         end
 570     end
 571 end
 572
 573 %%%% Make sure that each channel's intensity range is [0,1]
 574 maxV = sum(A.*(A>0),2);
 575 minV = sum(A.*(A<0),2);
 576 yNormal = (reshape(xRGB,[size(xRGB,1)*size(xRGB,2) 3]) * A' - repmat(minV, [1 size(xRGB,1)*size(xRGB,2)])') * diag(1./(maxV-minV)); % put in range [0,1]
 577 y = reshape(yNormal, [size(xRGB,1) size(xRGB,2) 3]);
 578
 579 %%%% The l2-norm of each of the 3 transform basis elements
 580 l2normLumChrom = diag(1./(maxV-minV))*sqrt(sum(A.^2,2));
 581
 582 if change_output
 583     A=l2normLumChrom;
 584 end
 585
 586 return;
 587
 588
 589
 590
 591 function yRGB=function_LumChrom2rgb(x,colormode)
 592 % Inverse color-space transformation   ( forward transformation is function_rgb2LumChrom.m )
 593 %
 594 % Alessandro Foi - Tampere University of Technology - 2005 - 2006   Public release v1.03 (March 2006)
 595 % -----------------------------------------------------------------------------------------------------------------------------------------------
 596 %
 597 % SYNTAX:
 598 %
 599 %   yRGB = function_LumChrom2rgb(x,colormode);
 600 %
 601 % INPUTS:
 602 %  x  is color-transformed image (with range typically included in or equal to [0 1]^3, depending on the transformation matrix)
 603 %
 604 %  colormode = 'opp', 'yCbCr', or a custom 3x3 matrix (e.g. provided by the forward transform when 'pca' is selected)
 605 %
 606 %       'opp'      opponent color space ('opp' is equirange version)
 607 %       'yCbCr'    standard yCbCr (e.g. for JPEG images)
 608 %
 609 % OUTPUTS:
 610 %   x  is RGB image (with range [0 1]^3)
 611 %
 612 %
 613 % NOTE:    'opp' is used by default if no colormode is specified
 614 %
 615
 616 if nargin==1
 617     colormode='opp';
 618 end
 619 if size(colormode)==[3 3]
 620     A=colormode;
 621     B=inv(A);
 622 else
 623     if strcmp(colormode,'opp')
 624         A =[1/3 1/3 1/3; 0.5  0  -0.5; 0.25  -0.5  0.25];
 625         B =[1 1 2/3;1 0 -4/3;1 -1 2/3];
 626     end
 627     if strcmp(colormode,'yCbCr')
 628         A=[0.299   0.587   0.114;   -0.16873660714285  -0.33126339285715   0.5;   0.5  -0.4186875  -0.0813125];
 629         B=inv(A);
 630     end
 631 end
 632
 633 %%%% Make sure that each channel's intensity range is [0,1]
 634 maxV = sum(A.*(A>0),2);
 635 minV = sum(A.*(A<0),2);
 636 xNormal = reshape(x,[size(x,1)*size(x,2) 3]) * diag(maxV-minV) +  repmat(minV, [1 size(x,1)*size(x,2)])'; % put in range [0,1]
 637 yRGB = reshape(xNormal * B', [ size(x,1) size(x,2) 3]);
 638
 639 return;
 640