BM3D/IDDBM3D/Demo_IDDBM3D.m

   1 function  [isnr, y_hat] = Demo_IDDBM3D(experiment_number, test_image_name)\r
   2 % ------------------------------------------------------------------------------------------\r
   3 %\r
   4 %     Demo software for BM3D-frame based image deblurring\r
   5 %               Public release ver. 0.8 (beta) (June 03, 2011)\r
   6 %\r
   7 % ------------------------------------------------------------------------------------------\r
   8 %\r
   9 %  This function implements the IDDBM3D image deblurring algorithm proposed in:\r
  10 %\r
  11 %  [1] A.Danielyan, V. Katkovnik, and K. Egiazarian, "BM3D frames and \r
  12 %   variational image deblurring," submitted to IEEE TIP, May 2011 \r
  13 %\r
  14 % ------------------------------------------------------------------------------------------\r
  15 %\r
  16 % authors:               Aram Danielyan\r
  17 %                        Vladimir Katkovnik\r
  18 %\r
  19 % web page:              http://www.cs.tut.fi/~foi/GCF-BM3D/\r
  20 %\r
  21 % contact:               firstname.lastname@tut.fi\r
  22 %\r
  23 % ------------------------------------------------------------------------------------------\r
  24 % Copyright (c) 2011 Tampere University of Technology.\r
  25 % All rights reserved.\r
  26 % This work should be used for nonprofit purposes only.\r
  27 % ------------------------------------------------------------------------------------------\r
  28 %\r
  29 % Disclaimer\r
  30 % ----------\r
  31 %\r
  32 % Any unauthorized use of these routines for industrial or profit-oriented activities is\r
  33 % expressively prohibited. By downloading and/or using any of these files, you implicitly\r
  34 % agree to all the terms of the TUT limited license (included in the file Legal_Notice.txt).\r
  35 % ------------------------------------------------------------------------------------------\r
  36 \r
  37 \r
  38 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  39 %  FUNCTION INTERFACE:\r
  40 %\r
  41 %  [psnr, y_hat] = Demo_IDDBM3D(experiment_number, test_image_name)\r
  42 %  \r
  43 %  INPUT:\r
  44 %   1) experiment_number: 1 -> PSF 1, sigma^2 = 2\r
  45 %                         2 -> PSF 1, sigma^2 = 8\r
  46 %                         3 -> PSF 2, sigma^2 = 0.308\r
  47 %                         4 -> PSF 3, sigma^2 = 49\r
  48 %                         5 -> PSF 4, sigma^2 = 4\r
  49 %                         6 -> PSF 5, sigma^2 = 64\r
  50 %                         7-13 -> experiments 7-13 are not described in [1].\r
  51 %                         see this file for the blur and noise parameters.\r
  52 %   2) test_image_name:   a valid filename of a grayscale test image\r
  53 %\r
  54 %  OUTPUT:\r
  55 %   1) isnr           the output improvement in SNR, dB\r
  56 %   2) y_hat:         the restored image\r
  57 %\r
  58 %  ! The function can work without any of the input arguments, \r
  59 %   in which case, the internal default ones are used !\r
  60 %   \r
  61 %   To run this demo functions within the BM3D package should be accessible to Matlab \r
  62 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
  63 addpath('../')\r
  64 \r
  65 if ~exist('experiment_number','var'), experiment_number=3; end\r
  66 if ~exist('test_image_name','var'), test_image_name='Cameraman256.png'; end\r
  67 \r
  68 filename=test_image_name;\r
  69 \r
  70 if 1 % \r
  71     initType = 'bm3ddeb'; %use output of the BM3DDEB to initialize the algorithm\r
  72 else\r
  73         initType = 'zeros'; %use zero image to initialize the algorithm\r
  74 end\r
  75 \r
  76 matchType = 'bm3ddeb'; %build groups using output of the BM3DDEB algorithm\r
  77 numIt = 200;\r
  78 \r
  79 fprintf('Experiment number: %d\n', experiment_number);\r
  80 fprintf('Image: %s\n', filename);\r
  81 \r
  82 %% ------- Generating bservation ---------------------------------------------\r
  83 disp('--- Generating observation ----');\r
  84 y=im2double(imread(filename));\r
  85 \r
  86 [yN,xN]=size(y);\r
  87 \r
  88 switch experiment_number\r
  89     case 1\r
  90         sigma=sqrt(2)/255; \r
  91         for x1=-7:7; for x2=-7:7; h(x1+8,x2+8)=1/(x1^2+x2^2+1); end, end; h=h./sum(h(:));\r
  92     case 2\r
  93         sigma=sqrt(8)/255;\r
  94         s1=0; for a1=-7:7; s1=s1+1; s2=0; for a2=-7:7; s2=s2+1; h(s1,s2)=1/(a1^2+a2^2+1); end, end;  h=h./sum(h(:));\r
  95     case 3 \r
  96         BSNR=40;\r
  97         sigma=-1; % if "sigma=-1", then the value of sigma depends on the BSNR\r
  98         h=ones(9); h=h./sum(h(:));\r
  99     case 4\r
 100         sigma=7/255;\r
 101         h=[1 4 6 4 1]'*[1 4 6 4 1]; h=h./sum(h(:));  % PSF\r
 102     case 5\r
 103         sigma=2/255;\r
 104         h=fspecial('gaussian', 25, 1.6);\r
 105     case 6\r
 106         sigma=8/255;\r
 107         h=fspecial('gaussian', 25, .4);\r
 108     %extra experiments\r
 109     case 7\r
 110         BSNR=30;\r
 111         sigma=-1;\r
 112         h=ones(9); h=h./sum(h(:));            \r
 113     case 8\r
 114         BSNR=20;\r
 115         sigma=-1;\r
 116         h=ones(9); h=h./sum(h(:));  \r
 117     case 9\r
 118         BSNR=40;\r
 119         sigma=-1;\r
 120         h=fspecial('gaussian', 25, 1.6);    \r
 121     case 10\r
 122         BSNR=20;\r
 123         sigma=-1;\r
 124         h=fspecial('gaussian', 25, 1.6);            \r
 125     case 11\r
 126         BSNR=15;\r
 127         sigma=-1; \r
 128         h=fspecial('gaussian', 25, 1.6);    \r
 129     case 12\r
 130         BSNR=40;\r
 131         sigma=-1; % if "sigma=-1", then the value of sigma depends on the BSNR\r
 132         h=ones(19); h=h./sum(h(:));            \r
 133     case 13\r
 134         BSNR=25;\r
 135         sigma=-1; % if "sigma=-1", then the value of sigma depends on the BSNR\r
 136         h=ones(19); h=h./sum(h(:));  \r
 137 end\r
 138 \r
 139 y_blur = imfilter(y, h, 'circular'); % performs blurring (by circular convolution)\r
 140 \r
 141 if sigma == -1;   %% check whether to use BSNR in order to define value of sigma\r
 142     sigma=sqrt(norm(y_blur(:)-mean(y_blur(:)),2)^2 /(yN*xN*10^(BSNR/10)));\r
 143     %     Xv% compute sigma from the desired BSNR\r
 144 end\r
 145 \r
 146 %%%% Create a blurred and noisy observation\r
 147 randn('seed',0);\r
 148 z = y_blur + sigma*randn(yN, xN);\r
 149 \r
 150 bsnr=10*log10(norm(y_blur(:)-mean(y_blur(:)),2)^2 /sigma^2/yN/xN);\r
 151 psnr_z =PSNR(y,z,1,0);\r
 152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\r
 153 \r
 154 fprintf('Observation BSNR: %4.2f, PSNR: %4.2f\n', bsnr, psnr_z);\r
 155 \r
 156 %% ----- Computing initial estimate ---------------------\r
 157 disp('--- Computing initial estimate  ----');\r
 158 \r
 159 [dummy, y_hat_RI,y_hat_RWI,zRI] = BM3DDEB_init(experiment_number, y, z, h, sigma);\r
 160 \r
 161 switch lower(initType)\r
 162     case 'zeros'\r
 163         y_hat_init=zeros(size(z));\r
 164     case 'zri'\r
 165         y_hat_init=zRI;\r
 166     case 'ri'\r
 167         y_hat_init=y_hat_RI;\r
 168     case 'bm3ddeb'\r
 169         y_hat_init=y_hat_RWI;\r
 170 \r
 171 end\r
 172 \r
 173 switch lower(matchType)\r
 174     case 'z'\r
 175         match_im = z;\r
 176     case 'y'\r
 177         match_im = y;\r
 178     case 'zri'\r
 179         match_im = zRI;\r
 180     case 'ri'\r
 181         match_im = y_hat_RI;\r
 182     case 'bm3ddeb'\r
 183         match_im = y_hat_RWI;   \r
 184 end\r
 185 \r
 186 psnr_init = PSNR(y, y_hat_init,1,0);\r
 187 \r
 188 fprintf('Initialization method: %s\n', initType);\r
 189 fprintf('Initial estimate ISNR: %4.2f, PSNR: %4.2f\n', psnr_init-psnr_z, psnr_init);\r
 190 \r
 191 %% ------- Core algorithm ---------------------\r
 192 %------ Description of the parameters of the IDDBM3D function ----------\r
 193 %y - true image (use [] if true image is unavaliable)\r
 194 %z - observed\r
 195 %h - blurring PSF\r
 196 %y_hat_init - initial estimate y_0\r
 197 %match_im - image used to constuct groups and calculate weights g_r\r
 198 %sigma - standard deviation of the noise\r
 199 %threshType = 'h'; %use 's' for soft thresholding\r
 200 %numIt - number of iterations\r
 201 %gamma - regularization parameter see [1]\r
 202 %tau - regularization parameter see [1] (thresholding level)\r
 203 %xi - regularization parameter see [1], it is always set to 1 in this implementation\r
 204 %showFigure - set to True to display figure with current estimate\r
 205 %--------------------------------------------------------------------\r
 206 \r
 207 threshType = 'h';\r
 208 showFigure = true;\r
 209 \r
 210 switch threshType\r
 211     case {'s'}\r
 212         gamma_tau_xi_inits= [\r
 213             0.0004509 0.70 1;%1\r
 214             0.0006803 0.78 1;%2\r
 215             0.0003485 0.65 1;%3\r
 216             0.0005259 0.72 1;%4\r
 217             0.0005327 0.82 1;%5\r
 218             7.632e-05 0.25 1;%6\r
 219             0.0005818 0.81 1;%7\r
 220             0.001149  1.18 1;%8\r
 221             0.0004155 0.74 1;%9\r
 222             0.0005591 0.74 1;%10\r
 223             0.0007989 0.82 1;%11\r
 224             0.0006702 0.75 1;%12\r
 225             0.001931  1.83 1;%13 \r
 226         ];\r
 227     case {'h'}\r
 228         gamma_tau_xi_inits= [ \r
 229             0.00051   3.13 1;%1\r
 230             0.0006004 2.75 1;%2\r
 231             0.0004573 2.91 1;%3\r
 232             0.0005959 2.82 1;%4\r
 233             0.0006018 3.63 1;%5\r
 234             0.0001726 2.24 1;%6\r
 235             0.00062   2.98 1;%7\r
 236             0.001047  3.80 1;%8\r
 237             0.0005125 3.00 1;%9\r
 238             0.0005685 2.80 1;%10\r
 239             0.0005716 2.75 1;%11\r
 240             0.0005938 2.55 1;%12\r
 241             0.001602  4.16 1;%13\r
 242         ];\r
 243 end\r
 244 \r
 245 gamma = gamma_tau_xi_inits(experiment_number,1);\r
 246 tau   = gamma_tau_xi_inits(experiment_number,2)/255*2.7;\r
 247 xi    = gamma_tau_xi_inits(experiment_number,3);\r
 248 \r
 249 disp('-------- Start ----------');\r
 250 fprintf('Number of iterations to perform: %d\n', numIt);\r
 251 fprintf('Thresholding type: %s\n', threshType);\r
 252 \r
 253 y_hat = IDDBM3D(y, h, z, y_hat_init, match_im, sigma, threshType, numIt, gamma, tau, xi, showFigure);\r
 254 psnr = PSNR(y,y_hat,1,0);\r
 255 isnr = psnr-psnr_z;\r
 256 \r
 257 disp('-------- Results --------');\r
 258 fprintf('Final estimate ISNR: %4.2f, PSNR: %4.2f\n', isnr, psnr);\r
 259 return;\r
 260 \r
 261 end\r
 262 \r
 263 function PSNRdb = PSNR(x, y, maxval, borders)\r
 264     if ~exist('borders', 'var'), borders = 0; end\r
 265     if ~exist('maxval', 'var'), maxval = 255; end\r
 266     \r
 267     xx=borders+1:size(x,1)-borders;\r
 268     yy=borders+1:size(x,2)-borders;\r
 269             \r
 270     PSNRdb = zeros(1,size(x,3));\r
 271     for fr=1:size(x,3) \r
 272         err = x(xx,yy,fr) - y(xx,yy,fr);\r
 273         PSNRdb(fr) = 10 * log10((maxval^2)/mean2(err.^2));    \r
 274     end\r
 275 end