idamnjanovic@10: function data=generateImageDenoiseProblem(im, trainnum, blocksize, dictsize, sigma, gain, maxval, initdict);
idamnjanovic@10: %%% Generate Image Denoising Problem
idamnjanovic@10: %
idamnjanovic@21: %   Centre for Digital Music, Queen Mary, University of London.
idamnjanovic@21: %   This file copyright 2010 Ivan Damnjanovic.
idamnjanovic@21: %
idamnjanovic@21: %   This program is free software; you can redistribute it and/or
idamnjanovic@21: %   modify it under the terms of the GNU General Public License as
idamnjanovic@21: %   published by the Free Software Foundation; either version 2 of the
idamnjanovic@21: %   License, or (at your option) any later version.  See the file
idamnjanovic@21: %   COPYING included with this distribution for more information.
idamnjanovic@21: %   
idamnjanovic@10: %   generateImageDenoiseProblem is a part of the SMALLbox and generates
idamnjanovic@10: %   a problem that can be used for comparison of Dictionary Learning/Sparse
idamnjanovic@10: %   Representation techniques in image denoising scenario.
idamnjanovic@10: %   The function takes as an input:
idamnjanovic@10: %   -   im - image matrix (if not present function promts user for an
idamnjanovic@10: %            image file) ,
idamnjanovic@10: %   -   trainnum - number of training samples (default - 40000)
idamnjanovic@10: %   -   blocksize - block (patch) vertical/horizontal dimension (default 8),
idamnjanovic@10: %   -   dictsize - dictionary size (default - 256),
idamnjanovic@10: %   -   sigma - noise level (default - 20),
idamnjanovic@10: %   -   noise gain (default - 1.15),
idamnjanovic@10: %   -   maxval - maximum value (default - 255)
idamnjanovic@10: %   -   initdict - initial dictionary (default - 4x overcomlete dct)
idamnjanovic@10: %
idamnjanovic@10: %   The output of the function is stucture with following fields:
idamnjanovic@10: %   -   name - name of the original image (if image is read inside of the
idamnjanovic@10: %              function)
idamnjanovic@10: %   -   Original - original image matrix,
idamnjanovic@10: %   -   Noisy - image with added noise,
idamnjanovic@10: %	-   b - training patches,
idamnjanovic@10: %	-   m - size of training patches (default 64),
idamnjanovic@10: %   -   n - number of training patches,
idamnjanovic@10: %   -   p - number of dictionary elements to be learned,
idamnjanovic@10: %   -   blocksize - block size (default [8 8]),
idamnjanovic@10: %   -   sigma - noise level,
idamnjanovic@10: %   -   noise gain (default - 1.15),
idamnjanovic@10: %   -   maxval - maximum value (default - 255)
idamnjanovic@10: %   -   initdict - initial dictionary (default - 4x overcomlete dct)
idamnjanovic@10: %   -   signalDim - signal dimension (default - 2)
idamnjanovic@10: %
idamnjanovic@10: %   Based on KSVD denoise demo by Ron Rubinstein
idamnjanovic@10: %   See also KSVDDENOISEDEMO and KSVDDEMO.
idamnjanovic@10: %   Ron Rubinstein
idamnjanovic@10: %   Computer Science Department
idamnjanovic@10: %   Technion, Haifa 32000 Israel
idamnjanovic@10: %   ronrubin@cs
idamnjanovic@10: %   August 2009
idamnjanovic@10: %%
idamnjanovic@10: disp(' ');
idamnjanovic@10: disp('  **********  Denoising Problem  **********');
idamnjanovic@10: disp(' ');
idamnjanovic@10: disp('  This function reads an image, adds random Gaussian noise,');
idamnjanovic@10: disp('  that can be later denoised by using dictionary learning techniques.');
idamnjanovic@10: disp(' ');
idamnjanovic@10: 
idamnjanovic@10: 
idamnjanovic@10: %% prompt user for image %%
idamnjanovic@10: %ask for file name
idamnjanovic@10: FS=filesep;
idamnjanovic@10: TMPpath=pwd;
idamnjanovic@10: if ~ exist( 'im', 'var' ) || isempty(im)
idamnjanovic@10:     [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
idamnjanovic@10:     cd([pathstr1,FS,'data',FS,'images']);
idamnjanovic@17:     [filename,pathname] = uigetfile({'*.png;'},'Select an image');
idamnjanovic@10:     [pathstr, name, ext, versn] = fileparts(filename);
idamnjanovic@10:     data.name=name;
idamnjanovic@10:     im = imread(filename);
idamnjanovic@44:     %im = double(im);
idamnjanovic@10: end;
idamnjanovic@44: im = double(im);
idamnjanovic@10: cd(TMPpath);
idamnjanovic@10: 
idamnjanovic@10: %% check input parameters %%
idamnjanovic@10: 
idamnjanovic@10: if ~ exist( 'blocksize', 'var' ) || isempty(blocksize),blocksize = 8;end
idamnjanovic@10: if ~ exist( 'dictsize', 'var' ) || isempty(dictsize), dictsize = 256;end
idamnjanovic@10: if ~ exist( 'trainnum', 'var' ) || isempty(trainnum),trainnum = 40000;end
idamnjanovic@10: if ~ exist( 'sigma', 'var' ) || isempty(sigma), sigma = 20; end
idamnjanovic@10: if ~ exist( 'gain', 'var' ) || isempty(gain), gain = 1.15; end
idamnjanovic@10: if ~ exist( 'maxval', 'var' ) || isempty(maxval), maxval = 255; end
idamnjanovic@10: if ~ exist( 'initdict', 'var' ) || isempty(initdict), initdict = 'odct'; end
idamnjanovic@10: 
idamnjanovic@10: %% generate noisy image %%
idamnjanovic@10: 
idamnjanovic@10: disp(' ');
idamnjanovic@10: disp('Generating noisy image...');
idamnjanovic@10: 
idamnjanovic@10: n = randn(size(im)) * sigma;
idamnjanovic@10: imnoise = im + n;
idamnjanovic@10: 
idamnjanovic@10: %% set parameters %%
idamnjanovic@10: 
idamnjanovic@10: x = imnoise;
idamnjanovic@10: p = ndims(x);
idamnjanovic@44: psnr=20*log10(maxval * sqrt(numel(im)) / norm(im(:)-imnoise(:)));
idamnjanovic@10: if (p==2 && any(size(x)==1) && length(blocksize)==1)
idamnjanovic@10:     p = 1;
idamnjanovic@10: end
idamnjanovic@10: 
idamnjanovic@10: % blocksize %
idamnjanovic@10: 
idamnjanovic@10: if (numel(blocksize)==1)
idamnjanovic@10:     blocksize = ones(1,p)*blocksize;
idamnjanovic@10: end
idamnjanovic@10: 
idamnjanovic@10: if (strcmpi(initdict,'odct'))
idamnjanovic@10:     initdict = odctndict(blocksize,dictsize,p);
idamnjanovic@10: elseif (strcmpi(initdict,'data'))
idamnjanovic@10:     clear initdict;    % causes initialization using random examples
idamnjanovic@10: else
idamnjanovic@10:     error('Invalid initial dictionary specified.');
idamnjanovic@10: end
idamnjanovic@10: 
idamnjanovic@10: if exist( 'initdict', 'var' )
idamnjanovic@10:     initdict = initdict(:,1:dictsize);
idamnjanovic@10: end
idamnjanovic@10: 
idamnjanovic@10: %%%% create training data %%%
idamnjanovic@10: 
idamnjanovic@10: ids = cell(p,1);
idamnjanovic@10: if (p==1)
idamnjanovic@10:     ids{1} = reggrid(length(x)-blocksize+1, trainnum, 'eqdist');
idamnjanovic@10: else
idamnjanovic@10:     [ids{:}] = reggrid(size(x)-blocksize+1, trainnum, 'eqdist');
idamnjanovic@10: end
idamnjanovic@10: X = sampgrid(x,blocksize,ids{:});
idamnjanovic@10: 
idamnjanovic@10: % remove dc in blocks to conserve memory %
idamnjanovic@10: 
idamnjanovic@10: bsize = 2000;
idamnjanovic@10: for i = 1:bsize:size(X,2)
idamnjanovic@10:     blockids = i : min(i+bsize-1,size(X,2));
idamnjanovic@10:     X(:,blockids) = remove_dc(X(:,blockids),'columns');
idamnjanovic@10: end
idamnjanovic@10: 
idamnjanovic@44: %  Noisy image blocks 
idamnjanovic@44: xcol=im2col(x,blocksize,'sliding');
idamnjanovic@44: [b1, dc] = remove_dc(xcol,'columns');
idamnjanovic@44: 
idamnjanovic@10: %% output structure %%
idamnjanovic@10: 
idamnjanovic@10: data.Original = im;
idamnjanovic@10: data.Noisy = imnoise;
idamnjanovic@44: data.noisy_psnr=psnr;
idamnjanovic@10: data.b = X;
idamnjanovic@44: data.b1=b1;
idamnjanovic@44: data.b1dc=dc;
idamnjanovic@10: data.m = size(X,1);
idamnjanovic@10: data.n = size(X,2);
idamnjanovic@10: data.p = dictsize;
idamnjanovic@10: data.blocksize=blocksize;
idamnjanovic@10: data.sigma = sigma;
idamnjanovic@10: data.gain = gain;
idamnjanovic@10: data.maxval = maxval;
idamnjanovic@10: data.initdict= initdict;
idamnjanovic@10: data.signalDim=2;
idamnjanovic@65: data.sparse=1;
idamnjanovic@10: end %% end of function