ivan@152: %% Dictionary Learning for Image Denoising - KSVD vs Recursive Least Squares ivan@152: % ivan@152: % This file contains an example of how SMALLbox can be used to test different ivan@152: % dictionary learning techniques in Image Denoising problem. ivan@152: % It calls generateImageDenoiseProblem that will let you to choose image, ivan@152: % add noise and use noisy image to generate training set for dictionary ivan@152: % learning. ivan@152: % Two dictionary learning techniques were compared: ivan@152: % - KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient ivan@152: % Implementation of the K-SVD Algorithm using Batch Orthogonal ivan@152: % Matching Pursuit", Technical Report - CS, Technion, April 2008. ivan@152: % - RLS-DLA - Skretting, K.; Engan, K.; , "Recursive Least Squares ivan@152: % Dictionary Learning Algorithm," Signal Processing, IEEE Transactions on, ivan@152: % vol.58, no.4, pp.2121-2130, April 2010 ivan@152: % ivan@152: ivan@152: ivan@152: % Centre for Digital Music, Queen Mary, University of London. ivan@152: % This file copyright 2011 Ivan Damnjanovic. ivan@152: % ivan@152: % This program is free software; you can redistribute it and/or ivan@152: % modify it under the terms of the GNU General Public License as ivan@152: % published by the Free Software Foundation; either version 2 of the ivan@152: % License, or (at your option) any later version. See the file ivan@152: % COPYING included with this distribution for more information. ivan@152: % ivan@152: %% ivan@152: ivan@152: ivan@152: ivan@152: % If you want to load the image outside of generateImageDenoiseProblem ivan@152: % function uncomment following lines. This can be useful if you want to ivan@152: % denoise more then one image for example. ivan@152: % Here we are loading test_image.mat that contains structure with 5 images : lena, ivan@152: % barbara,boat, house and peppers. ivan@152: clear; ivan@152: TMPpath=pwd; ivan@152: FS=filesep; ivan@152: [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); ivan@152: cd([pathstr1,FS,'data',FS,'images']); ivan@152: load('test_image.mat'); ivan@152: cd(TMPpath); ivan@152: ivan@152: % Deffining the noise levels that we want to test ivan@152: ivan@152: noise_level=[10 20 25 50 100]; ivan@152: ivan@152: % Here we loop through different noise levels and images ivan@152: ivan@153: for noise_ind=4:4 ivan@152: for im_num=1:1 ivan@152: ivan@152: % Defining Image Denoising Problem as Dictionary Learning ivan@152: % Problem. As an input we set the number of training patches. ivan@152: ivan@152: SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',256, noise_level(noise_ind)); ivan@152: SMALL.Problem.name=int2str(im_num); ivan@152: ivan@152: Edata=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain; ivan@152: maxatoms = floor(prod(SMALL.Problem.blocksize)/2); ivan@152: ivan@152: % results structure is to store all results ivan@152: ivan@152: results(noise_ind,im_num).noisy_psnr=SMALL.Problem.noisy_psnr; ivan@152: ivan@152: %% ivan@152: % Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary ivan@152: ivan@152: % Initialising Dictionary structure ivan@152: % Setting Dictionary structure fields (toolbox, name, param, D and time) ivan@152: % to zero values ivan@152: ivan@152: SMALL.DL(1)=SMALL_init_DL(); ivan@152: ivan@152: % Defining the parameters needed for dictionary learning ivan@152: ivan@152: SMALL.DL(1).toolbox = 'KSVD'; ivan@152: SMALL.DL(1).name = 'ksvd'; ivan@152: ivan@152: % Defining the parameters for KSVD ivan@152: % In this example we are learning 256 atoms in 20 iterations, so that ivan@152: % every patch in the training set can be represented with target error in ivan@152: % L2-norm (Edata) ivan@152: % Type help ksvd in MATLAB prompt for more options. ivan@152: ivan@152: ivan@152: SMALL.DL(1).param=struct(... ivan@152: 'Edata', Edata,... ivan@152: 'initdict', SMALL.Problem.initdict,... ivan@152: 'dictsize', SMALL.Problem.p,... ivan@152: 'exact', 1, ... ivan@152: 'iternum', 20,... ivan@152: 'memusage', 'high'); ivan@152: ivan@152: % Learn the dictionary ivan@152: ivan@152: SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1)); ivan@152: ivan@152: % Set SMALL.Problem.A dictionary ivan@152: % (backward compatiblity with SPARCO: solver structure communicate ivan@152: % only with Problem structure, ie no direct communication between DL and ivan@152: % solver structures) ivan@152: ivan@152: SMALL.Problem.A = SMALL.DL(1).D; ivan@161: SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); ivan@152: ivan@152: %% ivan@152: % Initialising solver structure ivan@152: % Setting solver structure fields (toolbox, name, param, solution, ivan@152: % reconstructed and time) to zero values ivan@152: ivan@152: SMALL.solver(1)=SMALL_init_solver; ivan@152: ivan@152: % Defining the parameters needed for image denoising ivan@152: ivan@152: SMALL.solver(1).toolbox='ompbox'; ivan@152: SMALL.solver(1).name='omp2'; ivan@152: SMALL.solver(1).param=struct(... ivan@152: 'epsilon',Edata,... ivan@152: 'maxatoms', maxatoms); ivan@152: ivan@152: % Denoising the image - find the sparse solution in the learned ivan@152: % dictionary for all patches in the image and the end it uses ivan@152: % reconstruction function to reconstruct the patches and put them into a ivan@152: % denoised image ivan@152: ivan@152: SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1)); ivan@152: ivan@152: % Show PSNR after reconstruction ivan@152: ivan@152: SMALL.solver(1).reconstructed.psnr ivan@152: ivan@152: %% ivan@152: % For comparison purposes we will denoise image with overcomplete DCT ivan@152: % here ivan@152: % Set SMALL.Problem.A dictionary to be oDCT (i.e. Problem.initdict - ivan@152: % since initial dictionaruy is already set to be oDCT when generating the ivan@152: % denoising problem ivan@152: ivan@152: ivan@152: % Initialising solver structure ivan@152: % Setting solver structure fields (toolbox, name, param, solution, ivan@152: % reconstructed and time) to zero values ivan@152: ivan@152: SMALL.solver(2)=SMALL_init_solver; ivan@152: ivan@152: % Defining the parameters needed for image denoising ivan@152: ivan@152: SMALL.solver(2).toolbox='ompbox'; ivan@152: SMALL.solver(2).name='omp2'; ivan@152: SMALL.solver(2).param=struct(... ivan@152: 'epsilon',Edata,... ivan@152: 'maxatoms', maxatoms); ivan@152: ivan@152: % Initialising Dictionary structure ivan@152: % Setting Dictionary structure fields (toolbox, name, param, D and time) ivan@152: % to zero values ivan@152: ivan@152: SMALL.DL(2)=SMALL_init_DL('TwoStepDL', 'MOD', '', 1); ivan@152: ivan@152: ivan@152: % Defining the parameters for MOD ivan@152: % In this example we are learning 256 atoms in 20 iterations, so that ivan@152: % every patch in the training set can be represented with target error in ivan@152: % L2-norm (EData) ivan@152: % Type help ksvd in MATLAB prompt for more options. ivan@152: ivan@152: ivan@152: SMALL.DL(2).param=struct(... ivan@152: 'solver', SMALL.solver(2),... ivan@152: 'initdict', SMALL.Problem.initdict,... ivan@152: 'dictsize', SMALL.Problem.p,... ivan@152: 'iternum', 40,... ivan@152: 'show_dict', 1); ivan@152: ivan@152: % Learn the dictionary ivan@152: ivan@152: SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2)); ivan@152: ivan@152: % Set SMALL.Problem.A dictionary ivan@152: % (backward compatiblity with SPARCO: solver structure communicate ivan@152: % only with Problem structure, ie no direct communication between DL and ivan@152: % solver structures) ivan@152: ivan@152: SMALL.Problem.A = SMALL.DL(2).D; ivan@161: SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); ivan@152: ivan@152: % Denoising the image - find the sparse solution in the learned ivan@152: % dictionary for all patches in the image and the end it uses ivan@152: % reconstruction function to reconstruct the patches and put them into a ivan@152: % denoised image ivan@152: ivan@152: SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2)); ivan@152: ivan@152: %% ivan@152: % In the b1 field all patches from the image are stored. For RLS-DLA we ivan@152: % will first exclude all the patches that have l2 norm smaller then ivan@152: % threshold and then take min(40000, number_of_remaining_patches) in ivan@152: % ascending order as our training set (SMALL.Problem.b) ivan@152: ivan@152: X=SMALL.Problem.b1; ivan@152: X_norm=sqrt(sum(X.^2, 1)); ivan@152: [X_norm_sort, p]=sort(X_norm); ivan@152: p1=p(X_norm_sort>Edata); ivan@152: if size(p1,2)>40000 ivan@152: p2 = randperm(size(p1,2)); ivan@152: p2=sort(p2(1:40000)); ivan@152: size(p2,2) ivan@152: SMALL.Problem.b=X(:,p1(p2)); ivan@152: else ivan@152: size(p1,2) ivan@152: SMALL.Problem.b=X(:,p1); ivan@152: ivan@152: end ivan@152: ivan@152: % Forgetting factor for RLS-DLA algorithm, in this case we are using ivan@152: % fixed value ivan@152: ivan@152: lambda=0.9998 ivan@152: ivan@152: % Use Recursive Least Squares ivan@152: % to Learn overcomplete dictionary ivan@152: ivan@152: % Initialising Dictionary structure ivan@152: % Setting Dictionary structure fields (toolbox, name, param, D and time) ivan@152: % to zero values ivan@152: ivan@152: SMALL.DL(3)=SMALL_init_DL(); ivan@152: ivan@152: % Defining fields needed for dictionary learning ivan@152: ivan@152: SMALL.DL(3).toolbox = 'SMALL'; ivan@152: SMALL.DL(3).name = 'SMALL_rlsdla'; ivan@152: SMALL.DL(3).param=struct(... ivan@152: 'Edata', Edata,... ivan@152: 'initdict', SMALL.Problem.initdict,... ivan@152: 'dictsize', SMALL.Problem.p,... ivan@152: 'forgettingMode', 'FIX',... ivan@152: 'forgettingFactor', lambda,... ivan@152: 'show_dict', 1000); ivan@152: ivan@152: ivan@152: SMALL.DL(3) = SMALL_learn(SMALL.Problem, SMALL.DL(3)); ivan@152: ivan@152: % Initialising solver structure ivan@152: % Setting solver structure fields (toolbox, name, param, solution, ivan@152: % reconstructed and time) to zero values ivan@152: ivan@152: SMALL.Problem.A = SMALL.DL(3).D; ivan@161: SMALL.Problem.reconstruct = @(x) ImageDenoise_reconstruct(x, SMALL.Problem); ivan@152: ivan@152: SMALL.solver(3)=SMALL_init_solver; ivan@152: ivan@152: % Defining the parameters needed for image denoising ivan@152: ivan@152: SMALL.solver(3).toolbox='ompbox'; ivan@152: SMALL.solver(3).name='omp2'; ivan@152: SMALL.solver(3).param=struct(... ivan@152: 'epsilon',Edata,... ivan@152: 'maxatoms', maxatoms); ivan@152: ivan@152: ivan@152: SMALL.solver(3)=SMALL_solve(SMALL.Problem, SMALL.solver(3)); ivan@152: ivan@152: SMALL.solver(3).reconstructed.psnr ivan@152: ivan@152: ivan@152: % show results % ivan@152: ivan@152: SMALL_ImgDeNoiseResult(SMALL); ivan@152: ivan@152: results(noise_ind,im_num).psnr.ksvd=SMALL.solver(1).reconstructed.psnr; ivan@152: results(noise_ind,im_num).psnr.odct=SMALL.solver(2).reconstructed.psnr; ivan@152: results(noise_ind,im_num).psnr.rlsdla=SMALL.solver(3).reconstructed.psnr; ivan@152: results(noise_ind,im_num).vmrse.ksvd=SMALL.solver(1).reconstructed.vmrse; ivan@152: results(noise_ind,im_num).vmrse.odct=SMALL.solver(2).reconstructed.vmrse; ivan@152: results(noise_ind,im_num).vmrse.rlsdla=SMALL.solver(3).reconstructed.vmrse; ivan@152: results(noise_ind,im_num).ssim.ksvd=SMALL.solver(1).reconstructed.ssim; ivan@152: results(noise_ind,im_num).ssim.odct=SMALL.solver(2).reconstructed.ssim; ivan@152: results(noise_ind,im_num).ssim.rlsdla=SMALL.solver(3).reconstructed.ssim; ivan@152: ivan@152: results(noise_ind,im_num).time.ksvd=SMALL.solver(1).time+SMALL.DL(1).time; ivan@152: results(noise_ind,im_num).time.rlsdla.time=SMALL.solver(3).time+SMALL.DL(3).time; ivan@152: clear SMALL; ivan@152: end ivan@152: end ivan@152: % save results.mat results