ivan@161: function reconstructed=ImageDenoise_reconstruct(y, Problem, SparseDict)
ivan@161: %%  Image Denoising Problem reconstruction function
ivan@161: %   
ivan@161: %   This reconstruction function is using sparse representation y 
ivan@161: %   in dictionary Problem.A to reconstruct the patches of the denoised
ivan@161: %   image.
ivan@161: 
ivan@161: %
ivan@161: %   Centre for Digital Music, Queen Mary, University of London.
ivan@161: %   This file copyright 2009 Ivan Damnjanovic.
ivan@161: %
ivan@161: %   This program is free software; you can redistribute it and/or
ivan@161: %   modify it under the terms of the GNU General Public License as
ivan@161: %   published by the Free Software Foundation; either version 2 of the
ivan@161: %   License, or (at your option) any later version.  See the file
ivan@161: %   COPYING included with this distribution for more information.
ivan@161: %%
ivan@161: 
ivan@161: 
ivan@161: % stepsize %
ivan@161: if (isfield(Problem,'stepsize'))
ivan@161:   stepsize = Problem.stepsize;
ivan@161:   if (numel(stepsize)==1)
ivan@161:     stepsize = ones(1,2)*stepsize;
ivan@161:   end
ivan@161: else
ivan@161:   stepsize = ones(1,2);
ivan@161: end
ivan@161: if (any(stepsize<1))
ivan@161:   error('Invalid step size.');
ivan@161: end
ivan@161: 
ivan@161: % lambda %
ivan@161: if (isfield(Problem,'lambda'))
ivan@161:   lambda = Problem.lambda;
ivan@161: else
ivan@161:   lambda = Problem.maxval/(10*Problem.sigma);
ivan@161: end
ivan@161: if exist('SparseDict','var')&&(SparseDict==1)
ivan@161:     if issparse(Problem.A)
ivan@161:         A = Problem.A;
ivan@161:       else
ivan@161:         A = sparse(Problem.A);
ivan@161:       end
ivan@161:     cl_samp=add_dc(dictsep(Problem.basedict,A,y), Problem.b1dc,'columns');
ivan@161: else
ivan@161:     cl_samp=add_dc(Problem.A*y, Problem.b1dc,'columns');
ivan@161: end
ivan@161: %   combine the patches into reconstructed image
ivan@161: cl_im=col2imstep(cl_samp, size(Problem.Noisy), Problem.blocksize);
ivan@161: 
ivan@161: cnt = countcover(size(Problem.Noisy),Problem.blocksize,stepsize);
ivan@161: 
ivan@161: im = (cl_im+lambda*Problem.Noisy)./(cnt + lambda);
ivan@161: % y(y~=0)=1;
ivan@161: % numD=sum(y,2);
ivan@161: % nnzy=sum(y,1);
ivan@161: % figure(200);plot(sort(numD));
ivan@161: % figure(201);plot(sort(nnzy));
ivan@161: [v.RMSErn, v.RMSEcd, v.rn_im, v.cd_im]=SMALL_vmrse_type2(Problem.Original, Problem.Noisy, im);
ivan@161: %% output structure image+psnr %%
ivan@161: reconstructed.Image=im;
ivan@161: reconstructed.psnr = 20*log10(Problem.maxval * sqrt(numel(Problem.Original(:))) / norm(Problem.Original(:)-im(:)));
ivan@161: reconstructed.vmrse=v;
ivan@161: reconstructed.ssim=SMALL_ssim_index(Problem.Original, im);
ivan@161: end