BM3D/BM3DSHARP.m

   1 function [y_hat] = BM3DSHARP(z, sigma, alpha_sharp, profile, print_to_screen)
   2 %
   3 %  Joint sharpening and denoising with BM3D. This is implementation of the
   4 %  BM3D-SH3D sharpening method that is developed in:
   5 %
   6 %  [1] K. Dabov, A. Foi, V. Katkovnik, and K. Egiazarian, "Joint image
   7 %   sharpening and denoising by 3D transform-domain collaborative filtering,"
   8 %   Proc. 2007 Int. TICSP Workshop Spectral Meth. Multirate Signal Process.,
   9 %   SMMSP 2007, Moscow, Russia, September 2007.
  10 %
  11 %  FUNCTION INTERFACE:
  12 %
  13 %  [ysharp] = BM3DSHARP(z, sigma, alpha_sharp, profile, print_to_screen)
  14 %
  15 %  The function can work without any of the input arguments, hence they are
  16 %  optional!
  17 %
  18 %  INPUTS (OPTIONAL):
  19 %
  20 %        1) z (matrix, size MxN)       : Input image (noisy and with poor contrast)
  21 %        2) sigma (double)             : Noise (IF ANY noise) standard deviation (signal assumed
  22 %                                          in the range [0, 255])
  23 %        3) alpha_sharp (double)       : Sharpening parameter (default: 1.5):
  24 %                                         (1,inf) -> sharpen
  25 %                                          1      -> no sharpening
  26 %                                         (0,1)   -> de-sharpen
  27 %        4) profile (char vector)      : 'lc' --> fast
  28 %                                        'np' --> normal (default)
  29 %        5) print_to_screen (boolean)  : 0 --> do not print output
  30 %                                          information (and do not plot figures)
  31 %                                        1 --> print figures (default)
  32 %
  33 %   OUTPUTS:
  34 %        1) ysharp (matrix, size MxN)  : Sharpened image (in the range [0,1])
  35 %
  36 %  BASIC USAGE EXAMPLES:
  37 %
  38 %     sigma = 10;
  39 %     z = im2double(imread('cameraman.tif'));
  40 %     z = z + (sigma/255)*randn(size(z));
  41 %     alpha_sharp = 1.3;
  42 %     [ysharp] = BM3DSHARP(z, sigma, alpha_sharp);
  43 %
  44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  45 %
  46 % Copyright © 2007 Tampere University of Technology. All rights reserved.
  47 % This work should only be used for nonprofit purposes.
  48 %
  49 % AUTHORS:
  50 %     Kostadin Dabov (2007), email: kostadin.dabov _at_ tut.fi
  51 %
  52 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  53
  54
  55 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  56 %%%% In case, an input image z is not provided, then use the filename
  57 %%%%  below to read an original image (might contain path also). Later,
  58 %%%%  artificial AWGN noise is added and this noisy image is processed
  59 %%%%  by the BM3D-SH3D.
  60 %%%%
  61 if (exist('image_name') ~= 1)
  62     image_name = [
  63     %
  64     %%%% Grayscale images
  65     %     'barco.png'
  66     %     'pentagon.tif'
  67          'Cameraman256.png'
  68     %     'boat.png'
  69     %     'Lena512.png'
  70     %     'house.png'
  71     %     'barbara.png'
  72 ];
  73 end
  74 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  75 %%%%  Quality/complexity trade-off profile selection
  76 %%%%
  77 %%%%  'np' --> Normal Profile (balanced quality)
  78 %%%%  'lc' --> Low Complexity Profile (fast, lower quality)
  79 %%%%
  80 %%%%  'high' --> High Profile (high quality, not documented in [1])
  81 %%%%
  82 if (exist('profile') ~= 1)
  83     profile         = 'np'; %% default profile
  84 end
  85
  86 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  87 %%%%  Specify the std. dev. of the corrupting noise
  88 %%%%
  89 if (exist('sigma') ~= 1),
  90     if (exist('z') ~= 1)
  91         sigma = 20; %% default standard deviation of the AWGN
  92     else
  93         fprintf('Please specify value for the s.t.d. "sigma"\n');
  94         y_hat = 0;
  95         return;
  96     end
  97 end
  98
  99 if (exist('alpha_sharp') ~= 1)
 100     alpha_sharp = 3/2;
 101 end
 102
 103 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 104 %%%% Following are the parameters for the Normal Profile.
 105 %%%%
 106
 107 %%%% Select transforms ('dct', 'dst', 'hadamard', or anything that is listed by 'help wfilters'):
 108 transform_2D_HT_name     = 'bior1.5'; %% transform used for the HT filt. of size N1 x N1
 109 transform_3rd_dim_name   = 'haar';    %% transform used in the 3-rd dim, the same for HT and Wiener filt.
 110
 111 %%%% Hard-thresholding (HT) parameters:
 112 N1                  = 8;   %% N1 x N1 is the block size used for the hard-thresholding (HT) filtering
 113 Nstep               = 3;   %% sliding step to process every next reference block
 114 N2                  = 16;  %% maximum number of similar blocks (maximum size of the 3rd dimension of a 3D array)
 115 Ns                  = 39;  %% length of the side of the search neighborhood for full-search block-matching (BM), must be odd
 116 tau_match           = 3000;%% threshold for the block-distance (d-distance)
 117 lambda_thr2D        = 0;   %% threshold parameter for the coarse initial denoising used in the d-distance measure
 118 lambda_thr3D        = 2.7; %% threshold parameter for the hard-thresholding in 3D transform domain
 119 beta                = 2.0; %% parameter of the 2D Kaiser window used in the reconstruction
 120
 121 %%%% Block-matching parameters:
 122 stepFS              = 1;  %% step that forces to switch to full-search BM, "1" implies always full-search
 123 smallLN             = 'not used in np'; %% if stepFS > 1, then this specifies the size of the small local search neighb.
 124 thrToIncStep        = 8;  %% used in the HT filtering to increase the sliding step in uniform regions
 125
 126 if strcmp(profile, 'lc') == 1,
 127
 128     Nstep               = 6;
 129     Ns                  = 25;
 130
 131     thrToIncStep        = 3;
 132     smallLN             = 3;
 133     stepFS              = 6*Nstep;
 134
 135 end
 136
 137 if (strcmp(profile, 'vn') == 1) | (sigma > 40),
 138
 139     transform_2D_HT_name = 'dct';
 140
 141     N1                  = 12;
 142     Nstep               = 4;
 143
 144     lambda_thr3D        = 2.8;
 145     lambda_thr2D        = 2.0;
 146     thrToIncStep        = 3;
 147     tau_match           = 5000;
 148
 149 end
 150
 151 decLevel = 0;        %% dec. levels of the dyadic wavelet 2D transform for blocks (0 means full decomposition, higher values decrease the dec. number)
 152 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
 153
 154 if strcmp(profile, 'high') == 1, %% this profile is not documented in [1]
 155
 156     decLevel     = 1;
 157     Nstep        = 2;
 158     lambda_thr3D = 2.5;
 159     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
 160     thr_mask = vMask * vMask';
 161     beta         = 2.5;
 162     beta_wiener  = 1.5;
 163
 164 end
 165
 166 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 167 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 168 %%%% Note: touch below this point only if you know what you are doing!
 169 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 170
 171 %%% Check whether to dump information to the screen or remain silent
 172 dump_output_information = 1;
 173 if (exist('print_to_screen') == 1) & (print_to_screen == 0),
 174     dump_output_information = 0;
 175 end
 176
 177 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 178 %%%% Create transform matrices, etc.
 179 %%%%
 180 [Tfor, Tinv]   = getTransfMatrix(N1, transform_2D_HT_name, decLevel);     %% get (normalized) forward and inverse transform matrices
 181
 182 if (strcmp(transform_3rd_dim_name, 'haar') == 1) | (strcmp(transform_3rd_dim_name(end-2:end), '1.1') == 1),
 183     %%% If Haar is used in the 3-rd dimension, then a fast internal transform is used, thus no need to generate transform
 184     %%% matrices.
 185     hadper_trans_single_den         = {};
 186     inverse_hadper_trans_single_den = {};
 187 else
 188     %%% Create transform matrices. The transforms are later applied by
 189     %%% matrix-vector multiplication for the 1D case.
 190     for hpow = 0:ceil(log2(max(N2,N2_wiener))),
 191         h = 2^hpow;
 192         [Tfor3rd, Tinv3rd]   = getTransfMatrix(h, transform_3rd_dim_name, 0);
 193         hadper_trans_single_den{h}         = single(Tfor3rd);
 194         inverse_hadper_trans_single_den{h} = single(Tinv3rd');
 195     end
 196 end
 197
 198 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 199 %%%% 2D Kaiser windows used in the aggregation of block-wise estimates
 200 %%%%
 201 if beta==2 & N1==8 % hardcode the window function so that the signal processing toolbox is not needed by default
 202     Wwin2D = [ 0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924;
 203         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 204         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 205         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 206         0.4325    0.6717    0.8644    0.9718    0.9718    0.8644    0.6717    0.4325;
 207         0.3846    0.5974    0.7688    0.8644    0.8644    0.7688    0.5974    0.3846;
 208         0.2989    0.4642    0.5974    0.6717    0.6717    0.5974    0.4642    0.2989;
 209         0.1924    0.2989    0.3846    0.4325    0.4325    0.3846    0.2989    0.1924];
 210 else
 211     Wwin2D = kaiser(N1, beta) * kaiser(N1, beta)'; % Kaiser window used in the aggregation of the HT part
 212 end
 213 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 214 %%%% If needed, read images, generate noise, or scale the images to the
 215 %%%% [0,1] interval
 216 %%%%
 217 if (exist('z') ~= 1)
 218     y        = im2double(imread(image_name));  %% read a noise-free image and put in intensity range [0,1]
 219     randn('seed', 0);                          %% generate seed
 220     z        = y + (sigma/255)*randn(size(y)); %% create a noisy image
 221 else  % external images
 222
 223     image_name = 'External image';
 224
 225     % convert z to double precision if needed
 226     z = double(z);
 227
 228     % if z's range is [0, 255], then convert to [0, 1]
 229     if (max(z(:)) > 10), % a naive check for intensity range
 230         z = z / 255;
 231     end
 232
 233 end
 234
 235 if (size(z,3) ~= 1),
 236     error('BM3D-SH3D accepts only grayscale 2D images.');
 237 end
 238
 239 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 240 %%% Print image information to the screen
 241 %%%%
 242 if dump_output_information == 1,
 243     fprintf('Image: %s (%dx%d), sigma: %.1f\n', image_name, size(z,1), size(z,2), sigma);
 244 end
 245
 246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 247 %%%% Apply the filtering MEX-subroutine
 248 %%%%
 249 tic;
 250 y_hat = bm3d_thr_sharpen_var(z, hadper_trans_single_den, Nstep, N1, N2, lambda_thr2D,...
 251          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, 1/alpha_sharp );
 252 estimate_elapsed_time = toc;
 253
 254 if dump_output_information == 1,
 255     fprintf('SHARPENING COMPLETED (total time: %.1f sec)\n', ...
 256         estimate_elapsed_time);
 257     imshow(z); figure, imshow(double(y_hat));
 258 end
 259
 260 return;
 261
 262
 263
 264 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 265 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 266 % Some auxiliary functions
 267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 268 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 269
 270
 271
 272 function [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 273 %
 274 % Create forward and inverse transform matrices, which allow for perfect
 275 % reconstruction. The forward transform matrix is normalized so that the
 276 % l2-norm of each basis element is 1.
 277 %
 278 % [Tforward, Tinverse] = getTransfMatrix (N, transform_type, dec_levels)
 279 %
 280 %  INPUTS:
 281 %
 282 %   N               --> Size of the transform (for wavelets, must be 2^K)
 283 %
 284 %   transform_type  --> 'dct', 'dst', 'hadamard', or anything that is
 285 %                       listed by 'help wfilters' (bi-orthogonal wavelets)
 286 %                       'DCrand' -- an orthonormal transform with a DC and all
 287 %                       the other basis elements of random nature
 288 %
 289 %   dec_levels      --> If a wavelet transform is generated, this is the
 290 %                       desired decomposition level. Must be in the
 291 %                       range [0, log2(N)-1], where "0" implies
 292 %                       full decomposition.
 293 %
 294 %  OUTPUTS:
 295 %
 296 %   Tforward        --> (N x N) Forward transform matrix
 297 %
 298 %   Tinverse        --> (N x N) Inverse transform matrix
 299 %
 300
 301 if exist('dec_levels') ~= 1,
 302     dec_levels = 0;
 303 end
 304
 305 if N == 1,
 306     Tforward = 1;
 307 elseif strcmp(transform_type, 'hadamard') == 1,
 308     Tforward    = hadamard(N);
 309 elseif (N == 8) & strcmp(transform_type, 'bior1.5')==1 % hardcoded transform so that the wavelet toolbox is not needed to generate it
 310     Tforward =  [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 311        0.219417649252501   0.449283757993216   0.449283757993216   0.219417649252501  -0.219417649252501  -0.449283757993216  -0.449283757993216  -0.219417649252501;
 312        0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846  -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284;
 313       -0.083506045090284   0.083506045090284  -0.083506045090284   0.083506045090284   0.569359398342846   0.402347308162278  -0.402347308162278  -0.569359398342846;
 314        0.707106781186547  -0.707106781186547                   0                   0                   0                   0                   0                   0;
 315                        0                   0   0.707106781186547  -0.707106781186547                   0                   0                   0                   0;
 316                        0                   0                   0                   0   0.707106781186547  -0.707106781186547                   0                   0;
 317                        0                   0                   0                   0                   0                   0   0.707106781186547  -0.707106781186547];
 318 elseif (N == 8) & strcmp(transform_type, 'dct')==1 % hardcoded transform so that the signal processing toolbox is not needed to generate it
 319     Tforward = [ 0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274   0.353553390593274;
 320        0.490392640201615   0.415734806151273   0.277785116509801   0.097545161008064  -0.097545161008064  -0.277785116509801  -0.415734806151273  -0.490392640201615;
 321        0.461939766255643   0.191341716182545  -0.191341716182545  -0.461939766255643  -0.461939766255643  -0.191341716182545   0.191341716182545   0.461939766255643;
 322        0.415734806151273  -0.097545161008064  -0.490392640201615  -0.277785116509801   0.277785116509801   0.490392640201615   0.097545161008064  -0.415734806151273;
 323        0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274   0.353553390593274  -0.353553390593274  -0.353553390593274   0.353553390593274;
 324        0.277785116509801  -0.490392640201615   0.097545161008064   0.415734806151273  -0.415734806151273  -0.097545161008064   0.490392640201615  -0.277785116509801;
 325        0.191341716182545  -0.461939766255643   0.461939766255643  -0.191341716182545  -0.191341716182545   0.461939766255643  -0.461939766255643   0.191341716182545;
 326        0.097545161008064  -0.277785116509801   0.415734806151273  -0.490392640201615   0.490392640201615  -0.415734806151273   0.277785116509801  -0.097545161008064];
 327 elseif (N == 8) & strcmp(transform_type, 'dst')==1 % hardcoded transform so that the PDE toolbox is not needed to generate it
 328     Tforward = [ 0.161229841765317   0.303012985114696   0.408248290463863   0.464242826880013   0.464242826880013   0.408248290463863   0.303012985114696   0.161229841765317;
 329        0.303012985114696   0.464242826880013   0.408248290463863   0.161229841765317  -0.161229841765317  -0.408248290463863  -0.464242826880013  -0.303012985114696;
 330        0.408248290463863   0.408248290463863                   0  -0.408248290463863  -0.408248290463863                   0   0.408248290463863   0.408248290463863;
 331        0.464242826880013   0.161229841765317  -0.408248290463863  -0.303012985114696   0.303012985114696   0.408248290463863  -0.161229841765317  -0.464242826880013;
 332        0.464242826880013  -0.161229841765317  -0.408248290463863   0.303012985114696   0.303012985114696  -0.408248290463863  -0.161229841765317   0.464242826880013;
 333        0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863                   0   0.408248290463863  -0.408248290463863;
 334        0.303012985114696  -0.464242826880013   0.408248290463863  -0.161229841765317  -0.161229841765317   0.408248290463863  -0.464242826880013   0.303012985114696;
 335        0.161229841765317  -0.303012985114696   0.408248290463863  -0.464242826880013   0.464242826880013  -0.408248290463863   0.303012985114696  -0.161229841765317];
 336 elseif strcmp(transform_type, 'dct') == 1,
 337     Tforward    = dct(eye(N));
 338 elseif strcmp(transform_type, 'dst') == 1,
 339     Tforward    = dst(eye(N));
 340 elseif strcmp(transform_type, 'DCrand') == 1,
 341     x = randn(N); x(1:end,1) = 1; [Q,R] = qr(x);
 342     if (Q(1) < 0),
 343         Q = -Q;
 344     end;
 345     Tforward = Q';
 346 else %% a wavelet decomposition supported by 'wavedec'
 347     %%% Set periodic boundary conditions, to preserve bi-orthogonality
 348     dwtmode('per','nodisp');
 349
 350     Tforward = zeros(N,N);
 351     for i = 1:N
 352         Tforward(:,i)=wavedec(circshift([1 zeros(1,N-1)],[dec_levels i-1]), log2(N), transform_type);  %% construct transform matrix
 353     end
 354 end
 355
 356 %%% Normalize the basis elements
 357 Tforward = (Tforward' * diag(sqrt(1./sum(Tforward.^2,2))))';
 358
 359 %%% Compute the inverse transform matrix
 360 Tinverse = inv(Tforward);
 361
 362 return;
 363