Mercurial > hg > smallbox

changeset 126:db5a7fe1a404 ivand_dev

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge from branch "sup_158_IMG_Processing_toolbox_"
author Ivan Damnjanovic lnx <ivan.damnjanovic@eecs.qmul.ac.uk>
date Wed, 25 May 2011 15:30:34 +0100
parents 5ded5e2e7d07 (current diff) 002ec1b2ceff (diff)
children 6f78b069e541 8e660fd14774
files util/ssim_index.m util/vmrse_type2.m
diffstat 13 files changed, 455 insertions(+), 288 deletions(-) [+]
line wrap: on
line diff
--- a/DL/RLS-DLA/SMALL_rlsdla.m	Wed May 25 13:33:47 2011 +0100
+++ b/DL/RLS-DLA/SMALL_rlsdla.m	Wed May 25 15:30:34 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
--- a/Problems/ImgDenoise_reconstruct.m	Wed May 25 13:33:47 2011 +0100
+++ b/Problems/ImgDenoise_reconstruct.m	Wed May 25 15:30:34 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
--- a/Problems/Pierre_Problem.m	Wed May 25 13:33:47 2011 +0100
+++ b/Problems/Pierre_Problem.m	Wed May 25 15:30:34 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 %%
 
--- a/Problems/generateImageDenoiseProblem.m	Wed May 25 13:33:47 2011 +0100
+++ b/Problems/generateImageDenoiseProblem.m	Wed May 25 15:30:34 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 %%
--- a/examples/Pierre Villars/Pierre_Villars_Example.m	Wed May 25 13:33:47 2011 +0100
+++ b/examples/Pierre Villars/Pierre_Villars_Example.m	Wed May 25 15:30:34 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
 
--- a/examples/SMALL_solver_test.m	Wed May 25 13:33:47 2011 +0100
+++ b/examples/SMALL_solver_test.m	Wed May 25 15:30:34 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
--- a/util/Pierre_reconstruct.m	Wed May 25 13:33:47 2011 +0100
+++ b/util/Pierre_reconstruct.m	Wed May 25 15:30:34 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;
--- a/util/SMALL_ImgDeNoiseResult.m	Wed May 25 13:33:47 2011 +0100
+++ b/util/SMALL_ImgDeNoiseResult.m	Wed May 25 15:30:34 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; 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/SMALL_showdict.m	Wed May 25 15:30:34 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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/SMALL_ssim_index.m	Wed May 25 15:30:34 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
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/util/SMALL_vmrse_type2.m	Wed May 25 15:30:34 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 = ((orig<recon_int)&(recon_int<=corr));
+    RN2 = ((orig>recon_int)&(recon_int>=corr));
+    CD1 = ((orig<recon_int)&(recon_int>corr));
+    CD2 = ((orig>recon_int)&(recon_int<corr));
+    
+    RMSErn = sqrt(sum(sum((RN1+RN2).*(orig-recon).^2)))/512;
+    RMSEcd = sqrt(sum(sum((CD1+CD2).*(orig-recon).^2)))/512;
+    rn_im=RN1+RN2;
+    cd_im=CD1+CD2;
+    
+end
+
--- a/util/ssim_index.m	Wed May 25 13:33:47 2011 +0100
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,194 +0,0 @@
-function [mssim, ssim_map] = 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
-%
-%========================================================================
-
-
-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);	%
-   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
--- a/util/vmrse_type2.m	Wed May 25 13:33:47 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 = ((orig<recon_int)&(recon_int<=corr));
-    RN2 = ((orig>recon_int)&(recon_int>=corr));
-    CD1 = ((orig<recon_int)&(recon_int>corr));
-    CD2 = ((orig>recon_int)&(recon_int<corr));
-    
-    RMSErn = sqrt(sum(sum((RN1+RN2).*(orig-recon).^2)))/512;
-    RMSEcd = sqrt(sum(sum((CD1+CD2).*(orig-recon).^2)))/512;
-    rn_im=RN1+RN2;
-    cd_im=CD1+CD2;
-    
-end
-