# HG changeset patch # User Ivan Damnjanovic lnx # Date 1306334077 -3600 # Node ID 6f78b069e541a3d95f4f692cd108b026ea762775 # Parent 9793ce00d72908b6ebdac581c4bc26656f3d2034# Parent db5a7fe1a4043ec338134816a04e38cadec1d3aa Merge from branch "ivand_dev" after issues 158, 159, 163 resolved diff -r 9793ce00d729 -r 6f78b069e541 DL/RLS-DLA/SMALL_rlsdla.m --- a/DL/RLS-DLA/SMALL_rlsdla.m Mon May 23 12:26:33 2011 +0100 +++ b/DL/RLS-DLA/SMALL_rlsdla.m Wed May 25 15:34:37 2011 +0100 @@ -133,8 +133,8 @@ end if (show_dictionary) - dictimg = showdict(D,[8 8],round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); - figure(2); imshow(imresize(dictimg,2,'nearest')); + dictimg = SMALL_showdict(D,[8 8],round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); + figure(2); imagesc(dictimg);colormap(gray);axis off; axis image; end % Forgetting factor @@ -193,9 +193,9 @@ C = C - (alfa * u)* u'; if (show_dictionary &&(mod(i,show_iter)==0)) - dictimg = showdict(D,[8 8],... + dictimg = SMALL_showdict(D,[8 8],... round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); - figure(2); imshow(imresize(dictimg,2,'nearest')); + figure(2); imagesc(dictimg);colormap(gray);axis off; axis image; pause(0.02); end end diff -r 9793ce00d729 -r 6f78b069e541 Problems/ImgDenoise_reconstruct.m --- a/Problems/ImgDenoise_reconstruct.m Mon May 23 12:26:33 2011 +0100 +++ b/Problems/ImgDenoise_reconstruct.m Wed May 25 15:34:37 2011 +0100 @@ -56,10 +56,10 @@ % nnzy=sum(y,1); % figure(200);plot(sort(numD)); % figure(201);plot(sort(nnzy)); -[v.RMSErn, v.RMSEcd, v.rn_im, v.cd_im]=vmrse_type2(Problem.Original, Problem.Noisy, im); +[v.RMSErn, v.RMSEcd, v.rn_im, v.cd_im]=SMALL_vmrse_type2(Problem.Original, Problem.Noisy, im); %% output structure image+psnr %% reconstructed.Image=im; reconstructed.psnr = 20*log10(Problem.maxval * sqrt(numel(Problem.Original(:))) / norm(Problem.Original(:)-im(:))); reconstructed.vmrse=v; -reconstructed.ssim=ssim_index(Problem.Original, im); +reconstructed.ssim=SMALL_ssim_index(Problem.Original, im); end \ No newline at end of file diff -r 9793ce00d729 -r 6f78b069e541 Problems/Pierre_Problem.m --- a/Problems/Pierre_Problem.m Mon May 23 12:26:33 2011 +0100 +++ b/Problems/Pierre_Problem.m Wed May 25 15:34:37 2011 +0100 @@ -1,14 +1,5 @@ function data=Pierre_Problem(src, trg, blocksize, dictsize); -%%% Generate Pierre Problem -% -% Centre for Digital Music, Queen Mary, University of London. -% This file copyright 2010 Ivan Damnjanovic. -% -% This program is free software; you can redistribute it and/or -% modify it under the terms of the GNU General Public License as -% published by the Free Software Foundation; either version 2 of the -% License, or (at your option) any later version. See the file -% COPYING included with this distribution for more information. +%% Generate Pierre Problem % % Pierre_Problem is a part of the SMALLbox and generates the problem % suggested by Professor Pierre Vandergheynst on the SMALL meeting in @@ -38,6 +29,15 @@ % - sparse - if 1 SMALL_solve will keep solution matrix in sparse form, % due to memory constrains. +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2010 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. %% prompt user for images %% % ask for source file name @@ -91,7 +91,7 @@ %% %% create dictionary data %% -S=im2col(src,blocksize,'sliding'); +S=im2colstep(src,blocksize); for j= 1:size(S,2) S(:,j)=S(:,j)./norm(S(:,j)); @@ -99,7 +99,7 @@ %% create measurement matrix %% -T=im2col(trg,blocksize, 'distinct'); +T=im2colstep(trg,blocksize, blocksize); %% output structure %% diff -r 9793ce00d729 -r 6f78b069e541 Problems/generateImageDenoiseProblem.m --- a/Problems/generateImageDenoiseProblem.m Mon May 23 12:26:33 2011 +0100 +++ b/Problems/generateImageDenoiseProblem.m Wed May 25 15:34:37 2011 +0100 @@ -1,14 +1,5 @@ function data=generateImageDenoiseProblem(im, trainnum, blocksize, dictsize, sigma, gain, maxval, initdict); -%%% Generate Image Denoising Problem -% -% Centre for Digital Music, Queen Mary, University of London. -% This file copyright 2010 Ivan Damnjanovic. -% -% This program is free software; you can redistribute it and/or -% modify it under the terms of the GNU General Public License as -% published by the Free Software Foundation; either version 2 of the -% License, or (at your option) any later version. See the file -% COPYING included with this distribution for more information. +%% Generate Image Denoising Problem % % generateImageDenoiseProblem is a part of the SMALLbox and generates % a problem that can be used for comparison of Dictionary Learning/Sparse @@ -39,6 +30,16 @@ % - maxval - maximum value (default - 255) % - initdict - initial dictionary (default - 4x overcomlete dct) % - signalDim - signal dimension (default - 2) + +% +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2010 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. % % Based on KSVD denoise demo by Ron Rubinstein % See also KSVDDENOISEDEMO and KSVDDEMO. @@ -136,7 +137,7 @@ end % Noisy image blocks -xcol=im2col(x,blocksize,'sliding'); +xcol=im2colstep(x,blocksize); [b1, dc] = remove_dc(xcol,'columns'); %% output structure %% diff -r 9793ce00d729 -r 6f78b069e541 SMALLboxSetup.m --- a/SMALLboxSetup.m Mon May 23 12:26:33 2011 +0100 +++ b/SMALLboxSetup.m Wed May 25 15:34:37 2011 +0100 @@ -1,5 +1,5 @@ function SMALLboxSetup(varargin) -%%% SMALLboxSetup +%% SMALLboxSetup % % Will automatically download and install existing toolboxes % on sparse representations and dictionary learning diff -r 9793ce00d729 -r 6f78b069e541 examples/Pierre Villars/Pierre_Villars_Example.m --- a/examples/Pierre Villars/Pierre_Villars_Example.m Mon May 23 12:26:33 2011 +0100 +++ b/examples/Pierre Villars/Pierre_Villars_Example.m Wed May 25 15:34:37 2011 +0100 @@ -28,10 +28,10 @@ % Show original image and image that is used as a dictionary figure('Name', 'Original and Dictionary Image'); -subplot(1,2,1); imshow(SMALL.Problem.imageTrg/SMALL.Problem.maxval); -title('Original Image'); -subplot(1,2,2); imshow(SMALL.Problem.imageSrc/SMALL.Problem.maxval); -title('Dictionary image:'); +subplot(1,2,1); imagesc(SMALL.Problem.imageTrg/SMALL.Problem.maxval); +title('Original Image');colormap(gray);axis off; axis image; +subplot(1,2,2); imagesc(SMALL.Problem.imageSrc/SMALL.Problem.maxval); +title('Dictionary image:');colormap(gray);axis off; axis image; % Using ten different dictionary sizes. First dictionary will contain all % patches from the source image and last one will have only @@ -80,10 +80,10 @@ %% show reconstructed image %% figure('Name', sprintf('dictsize=%d', dictsize(1,i))); - imshow(SMALL.solver(i).reconstructed.image/SMALL.Problem.maxval); + imagesc(SMALL.solver(i).reconstructed.image/SMALL.Problem.maxval); title(sprintf('Reconstructed image, PSNR: %.2f dB in %.2f s',... SMALL.solver(i).reconstructed.psnr, SMALL.solver(i).time )); - + colormap(gray);axis off; axis image; end diff -r 9793ce00d729 -r 6f78b069e541 examples/SMALL_solver_test.m --- a/examples/SMALL_solver_test.m Mon May 23 12:26:33 2011 +0100 +++ b/examples/SMALL_solver_test.m Wed May 25 15:34:37 2011 +0100 @@ -1,6 +1,5 @@ function SMALL_solver_test -% Example test of solvers from different toolboxes on Sparco compressed -% sensing problems +%% Example test of solvers from different toolboxes on Sparco problem 6 % % The main purpose of this example is to show how to use SMALL structure % to solve SPARCO compressed sensing problems (1-11) and compare results diff -r 9793ce00d729 -r 6f78b069e541 util/Pierre_reconstruct.m --- a/util/Pierre_reconstruct.m Mon May 23 12:26:33 2011 +0100 +++ b/util/Pierre_reconstruct.m Wed May 25 15:34:37 2011 +0100 @@ -1,5 +1,11 @@ function reconstructed=Pierre_reconstruct(y, Problem) -%%% Pierre Villars Example - reconstruction function +%% Pierre Villars Example - reconstruction function +% +% using sparse representation y in dictionary Problem.A reconstruct the +% patches from the target image +% This example is based on the experiment suggested by Professor Pierre +% Vandergheynst on the SMALL meeting in Villars. + % % Centre for Digital Music, Queen Mary, University of London. % This file copyright 2009 Ivan Damnjanovic. @@ -10,17 +16,12 @@ % License, or (at your option) any later version. See the file % COPYING included with this distribution for more information. % -% This example is based on the experiment suggested by Professor Pierre -% Vandergheynst on the SMALL meeting in Villars. - -% using sparse representation y in dictionary Problem.A reconstruct the -% patches from the target image - +%% imout=Problem.A*y; % combine the patches into reconstructed image -im=col2im(imout,Problem.blocksize,size(Problem.imageTrg),'disctint'); +im=col2imstep(imout,size(Problem.imageTrg),Problem.blocksize,Problem.blocksize); % bound the pixel values to [0,255] range im(im<0)=0; diff -r 9793ce00d729 -r 6f78b069e541 util/SMALL_ImgDeNoiseResult.m --- a/util/SMALL_ImgDeNoiseResult.m Mon May 23 12:26:33 2011 +0100 +++ b/util/SMALL_ImgDeNoiseResult.m Wed May 25 15:34:37 2011 +0100 @@ -1,5 +1,11 @@ -function SMALL_ImgDeNoiseResult(SMALL) +function SMALL_ImgDeNoiseResult(SMALL) +%% Represents the results of Dictionary Learning for Image denoising % +% Function gets as input SMALL structure and plots Image Denoise +% results: Original Image, Noisy Image and for learned dictionaries and +% denoised images +% + % Centre for Digital Music, Queen Mary, University of London. % This file copyright 2010 Ivan Damnjanovic. % @@ -8,10 +14,7 @@ % published by the Free Software Foundation; either version 2 of the % License, or (at your option) any later version. See the file % COPYING included with this distribution for more information. -% -% Function gets as input SMALL structure and plots Image Denoise -% results: Original Image, Noisy Image and for learned dictionaries and -% denoised images +%% figure('Name', sprintf('Image %s (training set size- %d, sigma - %d)',SMALL.Problem.name, SMALL.Problem.n, SMALL.Problem.sigma)); @@ -37,7 +40,7 @@ else D = SMALL.DL(i-1).D; end - dictimg = showdict(D,SMALL.Problem.blocksize,... + dictimg = SMALL_showdict(D,SMALL.Problem.blocksize,... round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); subplot(2,m,m+i);imagesc(dictimg);axis off; axis image; diff -r 9793ce00d729 -r 6f78b069e541 util/SMALL_learn.m --- a/util/SMALL_learn.m Mon May 23 12:26:33 2011 +0100 +++ b/util/SMALL_learn.m Wed May 25 15:34:37 2011 +0100 @@ -1,5 +1,13 @@ function DL = SMALL_learn(Problem,DL) -%%% SMALL Dictionary Learning +%% SMALL Dictionary Learning +% +% Function gets as input Problem and Dictionary Learning (DL) structures +% In Problem structure field b with the training set needs to be defined +% In DL fields with name of the toolbox and solver, and parameters file +% for particular dictionary learning technique needs to be present. +% +% Outputs are Learned dictionary and time spent as a part of DL structure + % % Centre for Digital Music, Queen Mary, University of London. % This file copyright 2009 Ivan Damnjanovic. @@ -9,13 +17,6 @@ % published by the Free Software Foundation; either version 2 of the % License, or (at your option) any later version. See the file % COPYING included with this distribution for more information. -% -% Function gets as input Problem and Dictionary Learning (DL) structures -% In Problem structure field b with the training set needs to be defined -% In DL fields with name of the toolbox and solver, and parameters file -% for particular dictionary learning technique needs to be present. -% -% Outputs are Learned dictionary and time spent as a part of DL structure %% fprintf('\nStarting Dictionary Learning %s... \n', DL.name); @@ -56,19 +57,7 @@ for i = 1: size(D,2) D(:,i)=D(:,i)/norm(D(:,i)); end - elseif strcmpi(DL.toolbox,'mpv2') - X = Problem.b(:,1:1:40000); - jD0 = mpv2.SimpleMatrix(DL.param.D); - jDicLea = mpv2.DictionaryLearning(jD0, 1,1); -% !!!!! MAYBE lambda is not needed for ILS - NEED TO CHECK!!!! - jDicLea.setLambda('L', DL.param.lambda, DL.param.lambda, 1000); - jDicLea.setVerbose(2); - jDicLea.setORMP(16, 1e-6, DL.param.abs); - eval(['jDicLea.',DL.name,'( X(:), DL.param.iternum );']); - jD = jDicLea.getDictionary(); - - D = reshape(jD.getAll(), size(X,1), DL.param.K); - + % To introduce new dictionary learning technique put the files in % your Matlab path. Next, unique name for your toolbox needs % to be defined and also prefferd API for toolbox functions diff -r 9793ce00d729 -r 6f78b069e541 util/SMALL_showdict.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_showdict.m Wed May 25 15:34:37 2011 +0100 @@ -0,0 +1,119 @@ +function x = SMALL_showdict(D,sz,n,m,varargin) +%% SMALL_SHOWDICT Display a dictionary of image patches. +%% Reimplementation of showdict function from KSVD toolbox with Image +%% Processing toolbox dependecies removed +% +% SMALL_SHOWDICT(D,SZ,N,M) displays the contents of the dictionary D, whos +% columns are 2-D image patches (in column-major order). SZ = [SX SY] is +% the size of the image patches. SHOWDICT displays the atoms on an N x M +% grid. If there are more atoms in D then only the first N*M are +% displayed. +% +% SMALL_SHOWDICT(...,'lines') separates the dictionary atoms by black lines. +% SMALL_SHOWDICT(...,'whitelines') separates the dictionary atoms by white +% lines. +% +% SMALL_SHOWDICT(...,'linewidth',W) when used with either 'lines' or +% 'whitelines' sets the width of the lines to W pixels (default=1). +% +% SMALL_SHOWDICT(...,'highcontrast') increases the contrast of the figure by +% normalizing the intensity values of each atom individually to the range +% of [0,1] (the default behavior is to normalize the values of the entire +% figure to [0,1] as one image). Note that in this way, the relative +% intensities of the atoms are not maintained. +% +% X = SMALL_SHOWDICT(...) returns a bitmat of the dictionary image without +% displaying the figure. + + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +% +%% + +if (size(D,2) < n*m) + D = [D zeros(size(D,1),n*m-size(D,2))]; +end + + +%%% parse input arguments %%% + +linewidth = 1; +highcontrast = 0; +drawlines = 0; +linecolor = 0; + +for i = 1:length(varargin) + if (~ischar(varargin{i})) + continue; + end + switch(varargin{i}) + case 'highcontrast' + highcontrast = 1; + case 'lines' + drawlines = 1; + case 'whitelines' + drawlines = 1; + linecolor = 1; + case 'linewidth' + linewidth = varargin{i+1}; + end +end + + + +%%% create dictionary image %%% + + +if (drawlines) + + D = [D ; nan(sz(1)*linewidth,size(D,2))]; + sz(2) = sz(2)+linewidth; + x = col2imstep(D(:,1:n*m),[n m].*sz, sz, sz); + sz = [sz(2) sz(1)]; + D = im2colstep(x',sz, sz); + D = [D ; nan(sz(1)*linewidth,size(D,2))]; + sz(2) = sz(2)+linewidth; + x = col2imstep(D(:,1:n*m),[m n].*sz,sz,sz); + x = x'; + x = x(1:end-linewidth,1:end-linewidth); + + if (highcontrast) + for i = 0:n-1 + for j = 0:m-1 + x(i*sz(1)+1:i*sz(1)+sz(1)-linewidth, j*sz(2)+1:j*sz(2)+sz(2)-linewidth) = ... + imnormalize(x(i*sz(1)+1:i*sz(1)+sz(1)-linewidth, j*sz(2)+1:j*sz(2)+sz(2)-linewidth)); + end + end + else + x = imnormalize(x); + end + + x(isnan(x)) = linecolor; + +else + + x = col2imstep(D(:,1:n*m),[n m].*sz, sz, sz); + + if (highcontrast) + for i = 0:n-1 + for j = 0:m-1 + x(i*sz(1)+1:i*sz(1)+sz(1), j*sz(2)+1:j*sz(2)+sz(2)) = ... + imnormalize(x(i*sz(1)+1:i*sz(1)+sz(1), j*sz(2)+1:j*sz(2)+sz(2))); + end + end + else + x = imnormalize(x); + end +end + + +if (nargout==0) + imagesc(dictimg);colormap(gray);axis off; axis image; +end diff -r 9793ce00d729 -r 6f78b069e541 util/SMALL_ssim_index.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_ssim_index.m Wed May 25 15:34:37 2011 +0100 @@ -0,0 +1,238 @@ +function [mssim, ssim_map] = SMALL_ssim_index(img1, img2, K, window, L) + +%======================================================================== +%SSIM Index, Version 1.0 +%Copyright(c) 2003 Zhou Wang +%All Rights Reserved. +% +%The author is with Howard Hughes Medical Institute, and Laboratory +%for Computational Vision at Center for Neural Science and Courant +%Institute of Mathematical Sciences, New York University. +% +%---------------------------------------------------------------------- +%Permission to use, copy, or modify this software and its documentation +%for educational and research purposes only and without fee is hereby +%granted, provided that this copyright notice and the original authors' +%names appear on all copies and supporting documentation. This program +%shall not be used, rewritten, or adapted as the basis of a commercial +%software or hardware product without first obtaining permission of the +%authors. The authors make no representations about the suitability of +%this software for any purpose. It is provided "as is" without express +%or implied warranty. +%---------------------------------------------------------------------- +% +%This is an implementation of the algorithm for calculating the +%Structural SIMilarity (SSIM) index between two images. Please refer +%to the following paper: +% +%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image +%quality assessment: From error measurement to structural similarity" +%IEEE Transactios on Image Processing, vol. 13, no. 1, Jan. 2004. +% +%Kindly report any suggestions or corrections to zhouwang@ieee.org +% +%---------------------------------------------------------------------- +% +%Input : (1) img1: the first image being compared +% (2) img2: the second image being compared +% (3) K: constants in the SSIM index formula (see the above +% reference). defualt value: K = [0.01 0.03] +% (4) window: local window for statistics (see the above +% reference). default widnow is Gaussian given by +% window = fspecial('gaussian', 11, 1.5); +% (5) L: dynamic range of the images. default: L = 255 +% +%Output: (1) mssim: the mean SSIM index value between 2 images. +% If one of the images being compared is regarded as +% perfect quality, then mssim can be considered as the +% quality measure of the other image. +% If img1 = img2, then mssim = 1. +% (2) ssim_map: the SSIM index map of the test image. The map +% has a smaller size than the input images. The actual size: +% size(img1) - size(window) + 1. +% +%Default Usage: +% Given 2 test images img1 and img2, whose dynamic range is 0-255 +% +% [mssim ssim_map] = ssim_index(img1, img2); +% +%Advanced Usage: +% User defined parameters. For example +% +% K = [0.05 0.05]; +% window = ones(8); +% L = 100; +% [mssim ssim_map] = ssim_index(img1, img2, K, window, L); +% +%See the results: +% +% mssim %Gives the mssim value +% imshow(max(0, ssim_map).^4) %Shows the SSIM index map +% + +% Removed dependency on Image Processing Toolbox - +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +%======================================================================== + + +if (nargin < 2 || nargin > 5) + ssim_index = -Inf; + ssim_map = -Inf; + return; +end + +if (size(img1) ~= size(img2)) + ssim_index = -Inf; + ssim_map = -Inf; + return; +end + +[M N] = size(img1); + +if (nargin == 2) + if ((M < 11) || (N < 11)) + ssim_index = -Inf; + ssim_map = -Inf; + return + end +% window = fspecial('gaussian', 11, 1.5); % +% +% Image Processing toolboxdependency +% - fspecial function is from IMP toolbox +% 2 - ways to do +% a) +% [wx,wy] = meshgrid(-5:5, -5:5); +% window = exp(-(wx.*wx + wy.*wy)/4.5); +% window = window/sum(window(:)); +% b) + window = [1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06;... + 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;... + 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;... + 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;... + 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;... + 0.000273561160085806,0.00202135875836257,0.00957662749024029,0.0290912256485504,0.0566619704916846,0.0707622377639470,0.0566619704916846,0.0290912256485504,0.00957662749024029,0.00202135875836257,0.000273561160085806;... + 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;... + 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;... + 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;... + 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;... + 1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06]; + K(1) = 0.01; % default settings + K(2) = 0.03; % + L = 255; % +end + +if (nargin == 3) + if ((M < 11) || (N < 11)) + ssim_index = -Inf; + ssim_map = -Inf; + return + end +% window = fspecial('gaussian', 11, 1.5); % +% +% Image Processing toolboxdependency +% - fspecial function is from IMP toolbox +% 2 - ways to do +% a) +% [wx,wy] = meshgrid(-5:5, -5:5); +% window = exp(-(wx.*wx + wy.*wy)/4.5); +% window = window/sum(window(:)); +% b) + window = [1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06;... + 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;... + 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;... + 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;... + 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;... + 0.000273561160085806,0.00202135875836257,0.00957662749024029,0.0290912256485504,0.0566619704916846,0.0707622377639470,0.0566619704916846,0.0290912256485504,0.00957662749024029,0.00202135875836257,0.000273561160085806;... + 0.000219050652866017,0.00161857756253439,0.00766836382523672,0.0232944324734871,0.0453713590956603,0.0566619704916846,0.0453713590956603,0.0232944324734871,0.00766836382523672,0.00161857756253439,0.000219050652866017;... + 0.000112464355116679,0.000831005429087199,0.00393706926284679,0.0119597604100370,0.0232944324734871,0.0290912256485504,0.0232944324734871,0.0119597604100370,0.00393706926284679,0.000831005429087199,0.000112464355116679;... + 3.70224770827489e-05,0.000273561160085806,0.00129605559384320,0.00393706926284679,0.00766836382523672,0.00957662749024029,0.00766836382523672,0.00393706926284679,0.00129605559384320,0.000273561160085806,3.70224770827489e-05;... + 7.81441153305360e-06,5.77411251978637e-05,0.000273561160085806,0.000831005429087199,0.00161857756253439,0.00202135875836257,0.00161857756253439,0.000831005429087199,0.000273561160085806,5.77411251978637e-05,7.81441153305360e-06;... + 1.05756559815326e-06,7.81441153305360e-06,3.70224770827489e-05,0.000112464355116679,0.000219050652866017,0.000273561160085806,0.000219050652866017,0.000112464355116679,3.70224770827489e-05,7.81441153305360e-06,1.05756559815326e-06]; + L = 255; + if (length(K) == 2) + if (K(1) < 0 || K(2) < 0) + ssim_index = -Inf; + ssim_map = -Inf; + return; + end + else + ssim_index = -Inf; + ssim_map = -Inf; + return; + end +end + +if (nargin == 4) + [H W] = size(window); + if ((H*W) < 4 || (H > M) || (W > N)) + ssim_index = -Inf; + ssim_map = -Inf; + return + end + L = 255; + if (length(K) == 2) + if (K(1) < 0 || K(2) < 0) + ssim_index = -Inf; + ssim_map = -Inf; + return; + end + else + ssim_index = -Inf; + ssim_map = -Inf; + return; + end +end + +if (nargin == 5) + [H W] = size(window); + if ((H*W) < 4 || (H > M) || (W > N)) + ssim_index = -Inf; + ssim_map = -Inf; + return + end + if (length(K) == 2) + if (K(1) < 0 || K(2) < 0) + ssim_index = -Inf; + ssim_map = -Inf; + return; + end + else + ssim_index = -Inf; + ssim_map = -Inf; + return; + end +end + +C1 = (K(1)*L)^2; +C2 = (K(2)*L)^2; +window = window/sum(sum(window)); +img1 = double(img1); +img2 = double(img2); + +mu1 = filter2(window, img1, 'valid'); +mu2 = filter2(window, img2, 'valid'); +mu1_sq = mu1.*mu1; +mu2_sq = mu2.*mu2; +mu1_mu2 = mu1.*mu2; +sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq; +sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq; +sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2; + +if (C1 > 0 & C2 > 0) + ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2)); +else + numerator1 = 2*mu1_mu2 + C1; + numerator2 = 2*sigma12 + C2; + denominator1 = mu1_sq + mu2_sq + C1; + denominator2 = sigma1_sq + sigma2_sq + C2; + ssim_map = ones(size(mu1)); + index = (denominator1.*denominator2 > 0); + ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index)); + index = (denominator1 ~= 0) & (denominator2 == 0); + ssim_map(index) = numerator1(index)./denominator1(index); +end + +mssim = sum(ssim_map(:), [], 'double')/numel(ssim_map); + +return \ No newline at end of file diff -r 9793ce00d729 -r 6f78b069e541 util/SMALL_vmrse_type2.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/util/SMALL_vmrse_type2.m Wed May 25 15:34:37 2011 +0100 @@ -0,0 +1,43 @@ +function [RMSErn, RMSEcd, rn_im, cd_im] = SMALL_vmrse_type2(orig, corr, recon) +%% Implementation of VectorRMSE type2 +% +% +% Input: +% - Original image +% - Corrupted image +% - Reconstructed Image +% +% Output: +% - RMSErn - RMSE from residual noise (noise not completely removed) +% - RMSEcd - RMSE from collateral distortion - excessive filtering +% - rn_im - image of residual noise +% - cd_im - image of collateral distortion +% +% F. Russo, "New Method for Performance Evaluation of Grayscale Image +% Denoising filters", IEEE Signal Processing Letters, vol. 17, no. 5, +% pp.417-420, May 2010 + +% Centre for Digital Music, Queen Mary, University of London. +% This file copyright 2011 Ivan Damnjanovic. +% +% This program is free software; you can redistribute it and/or +% modify it under the terms of the GNU General Public License as +% published by the Free Software Foundation; either version 2 of the +% License, or (at your option) any later version. See the file +% COPYING included with this distribution for more information. +%% + + recon_int = round(recon); + + RN1 = ((origrecon_int)&(recon_int>=corr)); + CD1 = ((origcorr)); + CD2 = ((orig>recon_int)&(recon_int 5) - ssim_index = -Inf; - ssim_map = -Inf; - return; -end - -if (size(img1) ~= size(img2)) - ssim_index = -Inf; - ssim_map = -Inf; - return; -end - -[M N] = size(img1); - -if (nargin == 2) - if ((M < 11) || (N < 11)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - window = fspecial('gaussian', 11, 1.5); % - K(1) = 0.01; % default settings - K(2) = 0.03; % - L = 255; % -end - -if (nargin == 3) - if ((M < 11) || (N < 11)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - window = fspecial('gaussian', 11, 1.5); - L = 255; - if (length(K) == 2) - if (K(1) < 0 || K(2) < 0) - ssim_index = -Inf; - ssim_map = -Inf; - return; - end - else - ssim_index = -Inf; - ssim_map = -Inf; - return; - end -end - -if (nargin == 4) - [H W] = size(window); - if ((H*W) < 4 || (H > M) || (W > N)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - L = 255; - if (length(K) == 2) - if (K(1) < 0 || K(2) < 0) - ssim_index = -Inf; - ssim_map = -Inf; - return; - end - else - ssim_index = -Inf; - ssim_map = -Inf; - return; - end -end - -if (nargin == 5) - [H W] = size(window); - if ((H*W) < 4 || (H > M) || (W > N)) - ssim_index = -Inf; - ssim_map = -Inf; - return - end - if (length(K) == 2) - if (K(1) < 0 || K(2) < 0) - ssim_index = -Inf; - ssim_map = -Inf; - return; - end - else - ssim_index = -Inf; - ssim_map = -Inf; - return; - end -end - -C1 = (K(1)*L)^2; -C2 = (K(2)*L)^2; -window = window/sum(sum(window)); -img1 = double(img1); -img2 = double(img2); - -mu1 = filter2(window, img1, 'valid'); -mu2 = filter2(window, img2, 'valid'); -mu1_sq = mu1.*mu1; -mu2_sq = mu2.*mu2; -mu1_mu2 = mu1.*mu2; -sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq; -sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq; -sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2; - -if (C1 > 0 & C2 > 0) - ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2)); -else - numerator1 = 2*mu1_mu2 + C1; - numerator2 = 2*sigma12 + C2; - denominator1 = mu1_sq + mu2_sq + C1; - denominator2 = sigma1_sq + sigma2_sq + C2; - ssim_map = ones(size(mu1)); - index = (denominator1.*denominator2 > 0); - ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index)); - index = (denominator1 ~= 0) & (denominator2 == 0); - ssim_map(index) = numerator1(index)./denominator1(index); -end - -mssim = mean2(ssim_map); - -return \ No newline at end of file diff -r 9793ce00d729 -r 6f78b069e541 util/vmrse_type2.m --- a/util/vmrse_type2.m Mon May 23 12:26:33 2011 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,43 +0,0 @@ -function [RMSErn, RMSEcd, rn_im, cd_im] = vmrse_type2(orig, corr, recon) - -%%% Implementation of VectorRMSE type2 -% -% Centre for Digital Music, Queen Mary, University of London. -% This file copyright 2011 Ivan Damnjanovic. -% -% This program is free software; you can redistribute it and/or -% modify it under the terms of the GNU General Public License as -% published by the Free Software Foundation; either version 2 of the -% License, or (at your option) any later version. See the file -% COPYING included with this distribution for more information. -% -% Input: -% - Original image -% - Corrupted image -% - Reconstructed Image -% -% Output: -% - RMSErn - RMSE from residual noise (noise not completely removed) -% - RMSEcd - RMSE from collateral distortion - excessive filtering -% - rn_im - image of residual noise -% - cd_im - image of collateral distortion -% -% F. Russo, "New Method for Performance Evaluation of Grayscale Image -% Denoising filters", IEEE Signal Processing Letters, vol. 17, no. 5, -% pp.417-420, May 2010 -%% - - recon_int = round(recon); - - RN1 = ((origrecon_int)&(recon_int>=corr)); - CD1 = ((origcorr)); - CD2 = ((orig>recon_int)&(recon_int