changeset 6:f72603404233

(none)
author idamnjanovic
date Mon, 22 Mar 2010 10:45:01 +0000
parents f44689e95ea4
children 0151f1ea080d
files examples/Automatic Music Transcription/SMALL_AMT_DL_test.m examples/Automatic Music Transcription/SMALL_AMT_KSVD_Err_test.m examples/Automatic Music Transcription/SMALL_AMT_KSVD_Sparsity_test.m examples/Automatic Music Transcription/SMALL_AMT_SPAMS_test.m examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsSPAMS.m examples/Image Denoising/SMALL_ImgDenoise_DL_test_SPAMS_lambda.m examples/Image Denoising/SMALL_ImgDenoise_DL_test_Training_size.m examples/Pierre Villars/Pierre_Villars_Example.m examples/private/add_dc.m examples/private/addtocols.c examples/private/addtocols.m examples/private/addtocols.mexw32 examples/private/collincomb.c examples/private/collincomb.m examples/private/collincomb.mexw32 examples/private/countcover.m examples/private/diag_ids.m examples/private/dictdist.m examples/private/imnormalize.m examples/private/iswhole.m examples/private/make.m examples/private/normcols.m examples/private/printf.m examples/private/reggrid.m examples/private/remove_dc.m examples/private/rowlincomb.c examples/private/rowlincomb.m examples/private/rowlincomb.mexw32 examples/private/sampgrid.m examples/private/secs2hms.m examples/private/spdiag.m examples/private/timerclear.m examples/private/timereta.m examples/private/timerinit.m
diffstat 34 files changed, 1242 insertions(+), 1309 deletions(-) [+]
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Automatic Music Transcription/SMALL_AMT_DL_test.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,211 @@
+%%  DICTIONARY LEARNING FOR AUTOMATIC MUSIC TRANSCRIPTION EXAMPLE 1
+%   This file contains an example of how SMALLbox can be used to test diferent
+%   dictionary learning techniques in Automatic Music Transcription problem.
+%   It calls generateAMT_Learning_Problem that will let you to choose midi,
+%   wave or mat file to be transcribe. If file is midi it will be first
+%   converted to wave and original midi file will be used for comparison with
+%   results of dictionary learning and reconstruction.
+%   The function will generarte the Problem structure that is used to learn
+%   Problem.p notes spectrograms from training set Problem.b using
+%   dictionary learning technique defined in DL structure.
+%
+%   Ivan Damnjanovic 2010
+%%
+
+clear;
+
+
+%   Defining Automatic Transcription of Piano tune as Dictionary Learning
+%   Problem
+
+SMALL.Problem = generateAMT_Learning_Problem();
+
+%%
+%   Use KSVD Dictionary Learning Algorithm to Learn 88 notes (defined in
+%   SMALL.Problem.p) using sparsity constrain only
+
+%   Initialising Dictionary structure
+%   Setting Dictionary structure fields (toolbox, name, param, D and time) 
+%   to zero values
+
+SMALL.DL(1)=SMALL_init_DL();
+
+%   Defining fields needed for dictionary learning
+
+SMALL.DL(1).toolbox = 'KSVD';
+SMALL.DL(1).name = 'ksvd';
+%   Defining the parameters for KSVD
+%   In this example we are learning 88 atoms in 100 iterations, so that
+%   every frame in the training set can be represented with maximum 3
+%   dictionary elements. Type help ksvd in MATLAB prompt for more options.
+
+SMALL.DL(1).param=struct(...
+    'Tdata', 3,...
+    'dictsize', SMALL.Problem.p,...
+    'iternum', 100);
+
+% Learn the dictionary
+
+SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1));
+
+%   Set SMALL.Problem.A dictionary and reconstruction function 
+%   (backward compatiblity with SPARCO: solver structure communicate
+%   only with Problem structure, ie no direct communication between DL and
+%   solver structures)
+
+SMALL.Problem.A = SMALL.DL(1).D;
+SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem);
+
+%%
+%   Initialising solver structure
+%   Setting solver structure fields (toolbox, name, param, solution, 
+%   reconstructed and time) to zero values
+%   As an example, SPAMS (Julien Mairal 2009) implementation of LARS
+%   algorithm is used for representation of training set in the learned
+%   dictionary.
+
+SMALL.solver(1)=SMALL_init_solver;
+
+%   Defining the parameters needed for sparse representation
+
+SMALL.solver(1).toolbox='SPAMS';
+SMALL.solver(1).name='mexLasso';
+
+%   Here we use mexLasso mode=2, with lambda=2, lambda2=0 and positivity
+%   constrain (type 'help mexLasso' for more information about modes):
+%   
+%   min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2
+
+SMALL.solver(1).param=struct(...
+    'lambda', 2,...
+    'pos', 1,...
+    'mode', 2);
+
+%   Call SMALL_soolve to represent the signal in the given dictionary. 
+%   As a final command SMALL_solve will call above defined reconstruction
+%   function to reconstruct the training set (Problem.b) in the learned 
+%   dictionary (Problem.A)
+
+SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1));
+
+%%
+%   Analysis of the result of automatic music transcription. If groundtruth
+%   exists, we can compare transcribed notes and original and get usual
+%   True Positives, False Positives and False Negatives measures.
+
+if ~isempty(SMALL.Problem.notesOriginal)
+    AMT_res(1) = AMT_analysis(SMALL.Problem, SMALL.solver(1));
+end
+
+
+%%
+%   Here we solve the same problem using non-negative sparse coding with 
+%   SPAMS online dictionary learning (Julien Mairal 2009)
+%
+
+%   Initialising Dictionary structure
+%   Setting Dictionary structure fields (toolbox, name, param, D and time) 
+%   to zero values
+
+SMALL.DL(2)=SMALL_init_DL();
+
+
+%   Defining fields needed for dictionary learning
+
+SMALL.DL(2).toolbox = 'SPAMS';
+SMALL.DL(2).name = 'mexTrainDL';
+
+%   Type 'help mexTrainDL in MATLAB prompt for explanation of parameters.
+
+SMALL.DL(2).param=struct(...
+    'K', SMALL.Problem.p,...
+    'lambda', 3,...
+    'iter', 300,...
+    'posAlpha', 1,...
+    'posD', 1,...
+    'whiten', 0,...
+    'mode', 2);
+
+%   Learn the dictionary
+
+SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2));
+
+%   Set SMALL.Problem.A dictionary and reconstruction function 
+%   (backward compatiblity with SPARCO: solver structure communicate
+%   only with Problem structure, ie no direct communication between DL and
+%   solver structures)
+
+SMALL.Problem.A = SMALL.DL(2).D;
+SMALL.Problem.reconstruct=@(x) SMALL_midiGenerate(x, SMALL.Problem);
+
+%%
+%   Initialising solver structure
+%   Setting solver structure fields (toolbox, name, param, solution, 
+%   reconstructed and time) to zero values
+%   As an example, SPAMS (Julien Mairal 2009) implementation of LARS
+%   algorithm is used for representation of training set in the learned
+%   dictionary.
+
+SMALL.solver(2)=SMALL_init_solver;
+
+%   Defining the parameters needed for sparse representation
+
+SMALL.solver(2).toolbox='SPAMS';
+SMALL.solver(2).name='mexLasso';
+
+%   Here we use mexLasso mode=2, with lambda=3, lambda2=0 and positivity
+%   constrain (type 'help mexLasso' for more information about modes):
+%   
+%   min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2
+
+SMALL.solver(2).param=struct('lambda', 3, 'pos', 1, 'mode', 2);
+
+%   Call SMALL_soolve to represent the signal in the given dictionary. 
+%   As a final command SMALL_solve will call above defined reconstruction
+%   function to reconstruct the training set (Problem.b) in the learned 
+%   dictionary (Problem.A)
+
+SMALL.solver(2)=SMALL_solve(SMALL.Problem, SMALL.solver(2));
+
+%%
+%   Analysis of the result of automatic music transcription. If groundtruth
+%   exists, we can compare transcribed notes and original and get usual
+%   True Positives, False Positives and False Negatives measures.
+
+if ~isempty(SMALL.Problem.notesOriginal)
+    AMT_res(2) = AMT_analysis(SMALL.Problem, SMALL.solver(2));
+end
+
+%%
+% Plot results and save midi files
+
+if ~isempty(SMALL.Problem.notesOriginal)
+    figAMT = SMALL_AMT_plot(SMALL, AMT_res);
+else
+    figAMT = figure('Name', 'Automatic Music Transcription KSVD vs SPAMS');
+    subplot(2,1,1); plot(SMALL.solver(1).reconstructed.notes(:,5), SMALL.solver(1).reconstructed.notes(:,3), 'kd ');
+            title (sprintf('%s dictionary in %.2f s', SMALL.DL(1).name, SMALL.DL(1).time));
+                xlabel('Time');
+                    ylabel('Note Number');
+    subplot(2,1,2); plot(SMALL.solver(2).reconstructed.notes(:,5), SMALL.solver(2).reconstructed.notes(:,3), 'b* ');
+            title (sprintf('%s dictionary in %.2f s', SMALL.DL(2).name, SMALL.DL(2).time));
+                xlabel('Time');
+                    ylabel('Note Number');
+end
+
+FS=filesep;
+
+[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+cd([pathstr1,FS,'results']);
+
+[filename,pathname] = uiputfile({' *.mid;' },'Save KSVD result midi');
+if filename~=0 writemidi(SMALL.solver(1).reconstructed.midi, [pathname,FS,filename]);end
+
+[filename,pathname] = uiputfile({' *.mid;' },'Save SPAMS result midi');
+if filename~=0 writemidi(SMALL.solver(2).reconstructed.midi, [pathname,FS,filename]);end
+
+[filename,pathname] = uiputfile({' *.fig;' },'Save KSVD vs SPAMS AMT figure');
+if filename~=0 saveas(figAMT, [pathname,FS,filename]);end
+
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Automatic Music Transcription/SMALL_AMT_KSVD_Err_test.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,127 @@
+%% DICTIONARY LEARNING FOR AUTOMATIC MUSIC TRANSCRIPTION
+%   This file contains an example of how SMALLbox can be used to test diferent
+%   dictionary learning techniques in Automatic Music Transcription problem.
+%   It calls generateAMT_Learning_Problem that will let you to choose midi,
+%   wave or mat file to be transcribe. If file is midi it will be first
+%   converted to wave and original midi file will be used for comparison with
+%   results of dictionary learning and reconstruction.
+%   The function will generarte the Problem structure that is used to learn
+%   Problem.p notes spectrograms from training set Problem.b using
+%   dictionary learning technique defined in DL structure.
+%
+%   Ivan Damnjanovic 2010
+%%
+
+clear;
+
+
+% Defining Automatic Transcription of Piano tune as Dictionary Learning
+% Problem
+
+SMALL.Problem = generateAMT_Learning_Problem();
+TPmax=0;
+for i=1:5
+    %%
+    %   Use KSVD Dictionary Learning Algorithm to Learn 88 notes (defined in
+    %   SMALL.Problem.p) using sparsity constrain only
+    
+    %   Initialising Dictionary structure
+    %   Setting Dictionary structure fields (toolbox, name, param, D and time)
+    %   to zero values
+    
+    SMALL.DL(i)=SMALL_init_DL(i);
+    
+    % Defining the parameters needed for dictionary learning
+    
+    SMALL.DL(i).toolbox = 'KSVD';
+    SMALL.DL(i).name = 'ksvd';
+    
+    %   Defining the parameters for KSVD
+    %   In this example we are learning 88 atoms in 100 iterations, so that
+    %   every frame in the training set can be represented with maximum 10
+    %   dictionary elements. However, our aim here is to show how individual
+    %   parameters can be ested in the AMT problem. We test five different
+    %   values for residual error (Edata) in KSVD algorithm.
+    %   Type help ksvd in MATLAB prompt for more options.
+    
+    Edata(i)=8+i*2;
+    SMALL.DL(i).param=struct(...
+        'Edata', Edata(i),...
+        'dictsize', SMALL.Problem.p,...
+        'iternum', 100,...
+        'maxatoms', 10);
+    
+    % Learn the dictionary
+    
+    SMALL.DL(i) = SMALL_learn(SMALL.Problem, SMALL.DL(i));
+    
+    %   Set SMALL.Problem.A dictionary and reconstruction function
+    %   (backward compatiblity with SPARCO: solver structure communicate
+    %   only with Problem structure, ie no direct communication between DL and
+    %   solver structures)
+    
+    SMALL.Problem.A = SMALL.DL(i).D;
+    SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem);
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    %   As an example, SPAMS (Julien Mairal 2009) implementation of LARS
+    %   algorithm is used for representation of training set in the learned
+    %   dictionary.
+    
+    SMALL.solver(1)=SMALL_init_solver;
+    
+    % Defining the parameters needed for sparse representation
+    
+    SMALL.solver(1).toolbox='SPAMS';
+    SMALL.solver(1).name='mexLasso';
+    SMALL.solver(1).param=struct('lambda', 2, 'pos', 1, 'mode', 2);
+    
+    %Represent Training set in the learned dictionary
+    
+    SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1));
+    
+    %%
+    %   Analysis of the result of automatic music transcription. If groundtruth
+    %   exists, we can compare transcribed notes and original and get usual
+    %   True Positives, False Positives and False Negatives measures.
+    
+    AMT_res(i) = AMT_analysis(SMALL.Problem, SMALL.solver(1));
+    
+    if AMT_res(i).TP>TPmax
+        TPmax=AMT_res(i).TP;
+        BLmidi=SMALL.solver(1).reconstructed.midi;
+        max=i;
+    end
+    
+end %end of for loop
+
+%%
+% Plot results and save midi files
+
+figAMTbest=SMALL_AMT_plot(SMALL, AMT_res(max));
+
+resFig=figure('Name', 'Automatic Music Transcription KSVD Error TEST');
+
+subplot (3,1,1); plot(Edata(:), [AMT_res(:).TP], 'ro-');
+title('True Positives vs Edata');
+
+subplot (3,1,2); plot(Edata(:), [AMT_res(:).FN], 'ro-');
+title('False Negatives vs Edata');
+
+subplot (3,1,3); plot(Edata(:), [AMT_res(:).FP], 'ro-');
+title('False Positives vs Edata');
+
+FS=filesep;
+[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+cd([pathstr1,FS,'results']);
+[filename,pathname] = uiputfile({' *.mid;' },'Save midi');
+if filename~=0 writemidi(BLmidi, [pathname,FS,filename]);end
+[filename,pathname] = uiputfile({' *.fig;' },'Save figure TP/FN/FP vs lambda');
+if filename~=0 saveas(resFig, [pathname,FS,filename]);end
+
+[filename,pathname] = uiputfile({' *.fig;' },'Save BEST AMT figure');
+if filename~=0 saveas(figAMTbest, [pathname,FS,filename]);end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Automatic Music Transcription/SMALL_AMT_KSVD_Sparsity_test.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,138 @@
+%%  DICTIONARY LEARNING FOR AUTOMATIC MUSIC TRANSCRIPTION EXAMPLE 1
+%   This file contains an example of how SMALLbox can be used to test diferent
+%   dictionary learning techniques in Automatic Music Transcription problem.
+%   It calls generateAMT_Learning_Problem that will let you to choose midi,
+%   wave or mat file to be transcribe. If file is midi it will be first
+%   converted to wave and original midi file will be used for comparison with
+%   results of dictionary learning and reconstruction.
+%   The function will generarte the Problem structure that is used to learn
+%   Problem.p notes spectrograms from training set Problem.b using
+%   dictionary learning technique defined in DL structure.
+%
+%   Ivan Damnjanovic 2010
+%%
+
+clear;
+
+
+% Defining Automatic Transcription of Piano tune as Dictionary Learning
+% Problem
+
+SMALL.Problem = generateAMT_Learning_Problem();
+
+TPmax=0;
+
+for i=1:10
+    
+    %%
+    %   Use KSVD Dictionary Learning Algorithm to Learn 88 notes (defined in
+    %   SMALL.Problem.p) using sparsity constrain only
+    
+    %   Initialising Dictionary structure
+    %   Setting Dictionary structure fields (toolbox, name, param, D and time)
+    %   to zero values
+    
+    SMALL.DL(i)=SMALL_init_DL(i);
+    
+    %   Defining fields needed for dictionary learning
+    
+    SMALL.DL(i).toolbox = 'KSVD';
+    SMALL.DL(i).name = 'ksvd';
+    
+    %   Defining the parameters for KSVD
+    %   In this example we are learning 88 atoms in 100 iterations.
+    %   our aim here is to show how individual parameters can be tested in
+    %   the AMT problem. We test ten different values for sparity (Tdata)
+    %   in KSVD algorithm.
+    %   Type help ksvd in MATLAB prompt for more options.
+    Tdata(i)=i;
+    SMALL.DL(i).param=struct('Tdata', Tdata(i), 'dictsize', SMALL.Problem.p, 'iternum', 100);
+    
+    % Learn the dictionary
+    
+    SMALL.DL(i) = SMALL_learn(SMALL.Problem, SMALL.DL(i));
+    
+    %   Set SMALL.Problem.A dictionary and reconstruction function
+    %   (backward compatiblity with SPARCO: solver structure communicate
+    %   only with Problem structure, ie no direct communication between DL and
+    %   solver structures)
+    
+    SMALL.Problem.A = SMALL.DL(i).D;
+    SMALL.Problem.reconstruct = @(x) SMALL_midiGenerate(x, SMALL.Problem);
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    %   As an example, SPAMS (Julien Mairal 2009) implementation of LARS
+    %   algorithm is used for representation of training set in the learned
+    %   dictionary.
+    
+    SMALL.solver(1)=SMALL_init_solver;
+    
+    % Defining the parameters needed for sparse representation
+    
+    SMALL.solver(1).toolbox='SPAMS';
+    SMALL.solver(1).name='mexLasso';
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    %   As an example, SPAMS (Julien Mairal 2009) implementation of LARS
+    %   algorithm is used for representation of training set in the learned
+    %   dictionary.
+    
+    SMALL.solver(1).param=struct(...
+        'lambda', 2,...
+        'pos', 1,...
+        'mode', 2);
+    
+    %   Call SMALL_soolve to represent the signal in the given dictionary.
+    %   As a final command SMALL_solve will call above defined reconstruction
+    %   function to reconstruct the training set (Problem.b) in the learned
+    %   dictionary (Problem.A)
+    
+    
+    SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1));
+    
+    %%
+    %   Analysis of the result of automatic music transcription. If groundtruth
+    %   exists, we can compare transcribed notes and original and get usual
+    %   True Positives, False Positives and False Negatives measures.
+    
+    AMT_res(i) = AMT_analysis(SMALL.Problem, SMALL.solver(1));
+    if AMT_res(i).TP>TPmax
+        TPmax=AMT_res(i).TP;
+        BLmidi=SMALL.solver(1).reconstructed.midi;
+        max=i;
+    end
+end % end of for loop
+
+%%
+% Plot results and save midi files
+
+figAMTbest=SMALL_AMT_plot(SMALL, AMT_res(max));
+
+resFig=figure('Name', 'Automatic Music Transcription KSVD Sparsity TEST');
+
+subplot (3,1,1); plot(Tdata(:), [AMT_res(:).TP], 'ro-');
+title('True Positives vs Tdata');
+
+subplot (3,1,2); plot(Tdata(:), [AMT_res(:).FN], 'ro-');
+title('False Negatives vs Tdata');
+
+subplot (3,1,3); plot(Tdata(:), [AMT_res(:).FP], 'ro-');
+title('False Positives vs Tdata');
+
+FS=filesep;
+[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+cd([pathstr1,FS,'results']);
+[filename,pathname] = uiputfile({' *.mid;' },'Save midi');
+if filename~=0 writemidi(BLmidi, [pathname,FS,filename]);end
+[filename,pathname] = uiputfile({' *.fig;' },'Save figure TP/FN/FP vs Tdata');
+if filename~=0 saveas(resFig, [pathname,FS,filename]);end
+
+[filename,pathname] = uiputfile({' *.fig;' },'Save BEST AMT figure');
+if filename~=0 saveas(figAMTbest, [pathname,FS,filename]);end
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Automatic Music Transcription/SMALL_AMT_SPAMS_test.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,140 @@
+%%  DICTIONARY LEARNING FOR AUTOMATIC MUSIC TRANSCRIPTION EXAMPLE 1
+%   This file contains an example of how SMALLbox can be used to test diferent
+%   dictionary learning techniques in Automatic Music Transcription problem.
+%   It calls generateAMT_Learning_Problem that will let you to choose midi,
+%   wave or mat file to be transcribe. If file is midi it will be first
+%   converted to wave and original midi file will be used for comparison with
+%   results of dictionary learning and reconstruction.
+%   The function will generarte the Problem structure that is used to learn
+%   Problem.p notes spectrograms from training set Problem.b using
+%   dictionary learning technique defined in DL structure.
+%
+%   Ivan Damnjanovic 2010
+%%
+
+clear;
+
+
+% Defining Automatic Transcription of Piano tune as Dictionary Learning
+% Problem
+
+SMALL.Problem = generateAMT_Learning_Problem();
+TPmax=0;
+%%
+for i=1:10
+    %%
+    %   Solving AMT problem using non-negative sparse coding with
+    %   SPAMS online dictionary learning (Julien Mairal 2009)
+    %
+    
+    %   Initialising Dictionary structure
+    %   Setting Dictionary structure fields (toolbox, name, param, D and time)
+    %   to zero values
+    
+    SMALL.DL(i)=SMALL_init_DL();
+    
+    %   Defining fields needed for dictionary learning
+    
+    SMALL.DL(i).toolbox = 'SPAMS';
+    SMALL.DL(i).name = 'mexTrainDL';
+    
+    %   We test SPAMS for ten different values of parameter lambda
+    %   Type 'help mexTrainDL in MATLAB prompt for explanation of parameters.
+    
+    lambda(i)=1.4+0.2*i;
+    
+    SMALL.DL(i).param=struct(...
+        'K', SMALL.Problem.p,...
+        'lambda', lambda(i),...
+        'iter', 300,...
+        'posAlpha', 1,...
+        'posD', 1,...
+        'whiten', 0,...
+        'mode', 2);
+    
+    %   Learn the dictionary
+    
+    SMALL.DL(i) = SMALL_learn(SMALL.Problem, SMALL.DL(i));
+    
+    %   Set SMALL.Problem.A dictionary and reconstruction function
+    %   (backward compatiblity with SPARCO: solver structure communicate
+    %   only with Problem structure, ie no direct communication between DL and
+    %   solver structures)
+    
+    SMALL.Problem.A = SMALL.DL(i).D;
+    SMALL.Problem.reconstruct=@(x) SMALL_midiGenerate(x, SMALL.Problem);
+    
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    %   As an example, SPAMS (Julien Mairal 2009) implementation of LARS
+    %   algorithm is used for representation of training set in the learned
+    %   dictionary.
+    
+    SMALL.solver(1)=SMALL_init_solver;
+    
+    %   Defining the parameters needed for sparse representation
+    
+    SMALL.solver(1).toolbox='SPAMS';
+    SMALL.solver(1).name='mexLasso';
+    
+    %   Here we use mexLasso mode=2, with lambda=3, lambda2=0 and positivity
+    %   constrain (type 'help mexLasso' for more information about modes):
+    %
+    %   min_{alpha_i} (1/2)||x_i-Dalpha_i||_2^2 + lambda||alpha_i||_1 + (1/2)lambda2||alpha_i||_2^2
+    
+    SMALL.solver(1).param=struct(...
+        'lambda', 3,...
+        'pos', 1,...
+        'mode', 2);
+    
+    %   Call SMALL_soolve to represent the signal in the given dictionary.
+    %   As a final command SMALL_solve will call above defined reconstruction
+    %   function to reconstruct the training set (Problem.b) in the learned
+    %   dictionary (Problem.A)
+    
+    SMALL.solver(1)=SMALL_solve(SMALL.Problem, SMALL.solver(1));
+    
+    %%
+    %   Analysis of the result of automatic music transcription. If groundtruth
+    %   exists, we can compare transcribed notes and original and get usual
+    %   True Positives, False Positives and False Negatives measures.
+    
+    AMT_res(i) = AMT_analysis(SMALL.Problem, SMALL.solver(1));
+    if AMT_res(i).TP>TPmax
+        TPmax=AMT_res(i).TP;
+        BLmidi=SMALL.solver(1).reconstructed.midi;
+        writemidi(SMALL.solver(1).reconstructed.midi, ['testL',i,'.mid']);
+        max=i;
+    end
+end %end of for loop
+%%
+% Plot results and save midi files
+
+figAMTbest=SMALL_AMT_plot(SMALL, AMT_res(max));
+
+resFig=figure('Name', 'Automatic Music Transcription SPAMS lambda TEST');
+
+subplot (3,1,1); plot(lambda(:), [AMT_res(:).TP], 'ro-');
+title('True Positives vs lambda');
+
+subplot (3,1,2); plot(lambda(:), [AMT_res(:).FN], 'ro-');
+title('False Negatives vs lambda');
+
+subplot (3,1,3); plot(lambda(:), [AMT_res(:).FP], 'ro-');
+title('False Positives vs lambda');
+
+FS=filesep;
+[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+cd([pathstr1,FS,'results']);
+[filename,pathname] = uiputfile({' *.mid;' },'Save midi');
+if filename~=0 writemidi(BLmidi, [pathname,FS,filename]);end
+[filename,pathname] = uiputfile({' *.fig;' },'Save figure TP/FN/FP vs lambda');
+if filename~=0 saveas(resFig, [pathname,FS,filename]);end
+
+[filename,pathname] = uiputfile({' *.fig;' },'Save BEST AMT figure');
+if filename~=0 saveas(figAMTbest, [pathname,FS,filename]);end
+
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsSPAMS.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,229 @@
+%% DICTIONARY LEARNING FOR IMAGE DENOISING
+%   This file contains an example of how SMALLbox can be used to test different
+%   dictionary learning techniques in Image Denoising problem.
+%   It calls generateImageDenoiseProblem that will let you to choose image,
+%   add noise and use noisy image to generate training set for dictionary
+%   learning.
+%   Three dictionary learning techniques were compared:
+%   -   KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient
+%              Implementation of the K-SVD Algorithm using Batch Orthogonal
+%              Matching Pursuit", Technical Report - CS, Technion, April 2008.
+%   -   KSVDS - R. Rubinstein, M. Zibulevsky, and M. Elad, "Learning Sparse
+%               Dictionaries for Sparse Signal Approximation", Technical
+%               Report - CS, Technion, June 2009.
+%   -   SPAMS - J. Mairal, F. Bach, J. Ponce and G. Sapiro. Online
+%               Dictionary Learning for Sparse Coding. International
+%               Conference on Machine Learning,Montreal, Canada, 2009
+%
+%
+% Ivan Damnjanovic 2010
+%%
+
+clear;
+
+%   If you want to load the image outside of generateImageDenoiseProblem
+%   function uncomment following lines. This can be useful if you want to
+%   denoise more then one image for example.
+
+% TMPpath=pwd;
+% FS=filesep;
+% [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+% cd([pathstr1,FS,'data',FS,'images']);
+% [filename,pathname] = uigetfile({'*.png;'},'Select a file containin pre-calculated notes');
+% [pathstr, name, ext, versn] = fileparts(filename);
+% test_image = imread(filename);
+% test_image = double(test_image);
+% cd(TMPpath);
+% SMALL.Problem.name=name;
+
+
+% Defining Image Denoising Problem as Dictionary Learning
+% Problem. As an input we set the number of training patches.
+
+SMALL.Problem = generateImageDenoiseProblem('', 40000);
+
+
+%%
+%   Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary
+
+%   Initialising Dictionary structure
+%   Setting Dictionary structure fields (toolbox, name, param, D and time)
+%   to zero values
+
+SMALL.DL(1)=SMALL_init_DL();
+
+% Defining the parameters needed for dictionary learning
+
+SMALL.DL(1).toolbox = 'KSVD';
+SMALL.DL(1).name = 'ksvd';
+
+%   Defining the parameters for KSVD
+%   In this example we are learning 256 atoms in 20 iterations, so that
+%   every patch in the training set can be represented with target error in
+%   L2-norm (EData)
+%   Type help ksvd in MATLAB prompt for more options.
+
+Edata=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain;
+SMALL.DL(1).param=struct(...
+    'Edata', Edata,...
+    'initdict', SMALL.Problem.initdict,...
+    'dictsize', SMALL.Problem.p,...
+    'iternum', 20,...
+    'memusage', 'high');
+
+%   Learn the dictionary
+
+SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1));
+
+%   Set SMALL.Problem.A dictionary
+%   (backward compatiblity with SPARCO: solver structure communicate
+%   only with Problem structure, ie no direct communication between DL and
+%   solver structures)
+
+SMALL.Problem.A = SMALL.DL(1).D;
+
+
+%%
+%   Initialising solver structure
+%   Setting solver structure fields (toolbox, name, param, solution,
+%   reconstructed and time) to zero values
+
+SMALL.solver(1)=SMALL_init_solver;
+
+% Defining the parameters needed for image denoising
+
+SMALL.solver(1).toolbox='ompbox';
+SMALL.solver(1).name='ompdenoise';
+
+%   Denoising the image - SMALL_denoise function is similar to SMALL_solve,
+%   but backward compatible with KSVD definition of denoising
+
+SMALL.solver(1)=SMALL_denoise(SMALL.Problem, SMALL.solver(1));
+
+%%
+% Use KSVDS Dictionary Learning Algorithm to denoise image
+
+%   Initialising solver structure
+%   Setting solver structure fields (toolbox, name, param, solution,
+%   reconstructed and time) to zero values
+
+SMALL.DL(2)=SMALL_init_DL();
+
+% Defining the parameters needed for dictionary learning
+
+SMALL.DL(2).toolbox = 'KSVDS';
+SMALL.DL(2).name = 'ksvds';
+
+%   Defining the parameters for KSVDS
+%   In this example we are learning 256 atoms in 20 iterations, so that
+%   every patch in the training set can be represented with target error in
+%   L2-norm (EDataS). We also impose "double sparsity" - dictionary itself
+%   has to be sparse in the given base dictionary (Tdict - number of
+%   nonzero elements per atom).
+%   Type help ksvds in MATLAB prompt for more options.
+
+EdataS=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain;
+SMALL.DL(2).param=struct(...
+    'Edata', EdataS, ...
+    'Tdict', 6,...
+    'stepsize', 1,...
+    'dictsize', SMALL.Problem.p,...
+    'iternum', 20,...
+    'memusage', 'high');
+SMALL.DL(2).param.initA = speye(SMALL.Problem.p);
+SMALL.DL(2).param.basedict{1} = odctdict(8,16);
+SMALL.DL(2).param.basedict{2} = odctdict(8,16);
+
+% Learn the dictionary
+
+SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2));
+
+%   Set SMALL.Problem.A dictionary and SMALL.Problem.basedictionary
+%   (backward compatiblity with SPARCO: solver structure communicate
+%   only with Problem structure, ie no direct communication between DL and
+%   solver structures)
+
+SMALL.Problem.A = SMALL.DL(2).D;
+SMALL.Problem.basedict{1} = SMALL.DL(2).param.basedict{1};
+SMALL.Problem.basedict{2} = SMALL.DL(2).param.basedict{2};
+
+%%
+%   Initialising solver structure
+%   Setting solver structure fields (toolbox, name, param, solution,
+%   reconstructed and time) to zero values
+
+SMALL.solver(2)=SMALL_init_solver;
+
+% Defining the parameters needed for image denoising
+
+SMALL.solver(2).toolbox='ompsbox';
+SMALL.solver(2).name='ompsdenoise';
+
+%   Denoising the image - SMALL_denoise function is similar to SMALL_solve,
+%   but backward compatible with KSVD definition of denoising
+%   Pay attention that since implicit base dictionary is used, denoising
+%   can be much faster then using explicit dictionary in KSVD example.
+
+SMALL.solver(2)=SMALL_denoise(SMALL.Problem, SMALL.solver(2));
+
+%%
+%   Use SPAMS Online Dictionary Learning Algorithm
+%   to Learn overcomplete dictionary (Julien Mairal 2009)
+%   (If you have not installed SPAMS please comment the following two cells)
+
+%   Initialising Dictionary structure
+%   Setting Dictionary structure fields (toolbox, name, param, D and time)
+%   to zero values
+
+SMALL.DL(3)=SMALL_init_DL();
+
+%   Defining fields needed for dictionary learning
+
+SMALL.DL(3).toolbox = 'SPAMS';
+SMALL.DL(3).name = 'mexTrainDL';
+
+%   Type 'help mexTrainDL in MATLAB prompt for explanation of parameters.
+
+SMALL.DL(3).param=struct(...
+    'D', SMALL.Problem.initdict,...
+    'K', SMALL.Problem.p,...
+    'lambda', 2,...
+    'iter', 200,...
+    'mode', 3, ...
+    'modeD', 0);
+
+%   Learn the dictionary
+
+SMALL.DL(3) = SMALL_learn(SMALL.Problem, SMALL.DL(3));
+
+%   Set SMALL.Problem.A dictionary
+%   (backward compatiblity with SPARCO: solver structure communicate
+%   only with Problem structure, ie no direct communication between DL and
+%   solver structures)
+
+SMALL.Problem.A = SMALL.DL(3).D;
+
+
+%%
+%   Initialising solver structure
+%   Setting solver structure fields (toolbox, name, param, solution,
+%   reconstructed and time) to zero values
+
+SMALL.solver(3)=SMALL_init_solver;
+
+% Defining the parameters needed for denoising
+
+SMALL.solver(3).toolbox='ompbox';
+SMALL.solver(3).name='ompdenoise';
+
+%   Denoising the image - SMALL_denoise function is similar to SMALL_solve,
+%   but backward compatible with KSVD definition of denoising
+
+SMALL.solver(3)=SMALL_denoise(SMALL.Problem, SMALL.solver(3));
+
+%%
+% Plot results and save midi files
+
+% show results %
+
+SMALL_ImgDeNoiseResult(SMALL);
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_SPAMS_lambda.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,121 @@
+%% DICTIONARY LEARNING FOR IMAGE DENOISING
+%   This file contains an example of how SMALLbox can be used to test different
+%   dictionary learning techniques in Image Denoising problem.
+%   This example can be used to test SPAMS for different values of
+%   parameter lambda. In no way it represents extensive testing of image
+%   denoising. It should only give an idea how SMALL structure can be used
+%   for testing.
+%
+%   Ivan Damnjanovic 2010
+%%
+
+clear all;
+
+%% Load an image
+TMPpath=pwd;
+FS=filesep;
+[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+cd([pathstr1,FS,'data',FS,'images']);
+[filename,pathname] = uigetfile({'*.png;'},'Select a file containin pre-calculated notes');
+[pathstr, name, ext, versn] = fileparts(filename);
+test_image = imread(filename);
+test_image = double(test_image);
+cd(TMPpath);
+%%
+
+%   number of different values we want to test
+
+n =4;
+
+lambda=zeros(1,n);
+time = zeros(2,n);
+psnr = zeros(2,n);
+
+for i=1:n
+    
+    %   Here we want to test time spent and quality of denoising for
+    %   different lambda parameters.
+    
+    lambda(i)=1+i*0.5;
+    
+    %   Defining Image Denoising Problem as Dictionary Learning Problem.
+    
+    SMALL.Problem = generateImageDenoiseProblem(test_image);
+    SMALL.Problem.name=name;
+    %%
+    %   Use SPAMS Online Dictionary Learning Algorithm
+    %   to Learn overcomplete dictionary (Julien Mairal 2009)
+    
+    %   Initialising Dictionary structure
+    %   Setting Dictionary structure fields (toolbox, name, param, D and time)
+    %   to zero values
+    
+    SMALL.DL(1)=SMALL_init_DL();
+    
+    %   Defining fields needed for dictionary learning
+    
+    SMALL.DL(1).toolbox = 'SPAMS';
+    SMALL.DL(1).name = 'mexTrainDL';
+    
+    %   Type 'help mexTrainDL in MATLAB prompt for explanation of parameters.
+    
+    SMALL.DL(1).param=struct(...
+        'D', SMALL.Problem.initdict,...
+        'K', SMALL.Problem.p,...
+        'lambda', lambda(i),...
+        'iter', 200,...
+        'mode', 3,...
+        'modeD', 0);
+    
+    %   Learn the dictionary
+    
+    SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1));
+    
+    %   Set SMALL.Problem.A dictionary
+    %   (backward compatiblity with SPARCO: solver structure communicate
+    %   only with Problem structure, ie no direct communication between DL and
+    %   solver structures)
+    
+    SMALL.Problem.A = SMALL.DL(1).D;
+    
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    
+    SMALL.solver(1)=SMALL_init_solver;
+    
+    % Defining the parameters needed for sparse representation
+    
+    SMALL.solver(1).toolbox='ompbox';
+    SMALL.solver(1).name='ompdenoise';
+    
+    %   Denoising the image - SMALL_denoise function is similar to SMALL_solve,
+    %   but backward compatible with KSVD definition of denoising
+    
+    SMALL.solver(1)=SMALL_denoise(SMALL.Problem, SMALL.solver(1));
+    
+    
+    %% show results %%
+    %   This will show denoised image and dictionary for all lambdas. If you
+    %   are not interested to see it and do not want clutter your screen
+    %   comment following line
+    
+    SMALL_ImgDeNoiseResult(SMALL);
+    
+    
+    time(1,i) = SMALL.DL(1).time;
+    psnr(1,i) = SMALL.solver(1).reconstructed.psnr;
+    
+    clear SMALL
+end
+
+%% show time and psnr %%
+figure('Name', 'SPAMS LAMBDA TEST');
+
+subplot(1,2,1); plot(lambda, time(1,:), 'ro-');
+title('time vs lambda');
+subplot(1,2,2); plot(lambda, psnr(1,:), 'b*-');
+title('PSNR vs lambda');
+
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_Training_size.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,189 @@
+%% DICTIONARY LEARNING FOR IMAGE DENOISING
+%   This file contains an example of how SMALLbox can be used to test different
+%   dictionary learning techniques in Image Denoising problem.
+%   It calls generateImageDenoiseProblem that will let you to choose image,
+%   add noise and use noisy image to generate training set for dictionary
+%   learning.
+%   We tested time and psnr for two dictionary learning techniques. This
+%   example does not represnt any extensive testing. The aim of this
+%   example is just to show how SMALL structure can be used for testing.
+%
+%   Two dictionary learning techniques were compared:
+%   -   KSVD - M. Elad, R. Rubinstein, and M. Zibulevsky, "Efficient
+%              Implementation of the K-SVD Algorithm using Batch Orthogonal
+%              Matching Pursuit", Technical Report - CS, Technion, April 2008.
+%   -   SPAMS - J. Mairal, F. Bach, J. Ponce and G. Sapiro. Online
+%               Dictionary Learning for Sparse Coding. International
+%               Conference on Machine Learning,Montreal, Canada, 2009
+%
+%
+% Ivan Damnjanovic 2010
+%%
+
+clear all;
+
+%% Load an image
+TMPpath=pwd;
+FS=filesep;
+[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m'));
+cd([pathstr1,FS,'data',FS,'images']);
+[filename,pathname] = uigetfile({'*.png;'},'Select a file containin pre-calculated notes');
+[pathstr, name, ext, versn] = fileparts(filename);
+test_image = imread(filename);
+test_image = double(test_image);
+cd(TMPpath);
+
+%   number of different values we want to test
+n =5;
+step = floor((size(test_image,1)-8+1)*(size(test_image,2)-8+1)/n);
+Training_size=zeros(1,n);
+time = zeros(2,n);
+psnr = zeros(2,n);
+for i=1:n
+    
+    %   Here we want to test time spent and quality of denoising for
+    %   different sizes of training sample.
+    Training_size(i)=i*step;
+    
+    SMALL.Problem = generateImageDenoiseProblem(test_image,Training_size(i));
+    SMALL.Problem.name=name;
+    %%
+    %   Use KSVD Dictionary Learning Algorithm to Learn overcomplete dictionary
+    
+    %   Initialising Dictionary structure
+    %   Setting Dictionary structure fields (toolbox, name, param, D and time)
+    %   to zero values
+    
+    SMALL.DL(1)=SMALL_init_DL();
+    
+    % Defining the parameters needed for dictionary learning
+    
+    SMALL.DL(1).toolbox = 'KSVD';
+    SMALL.DL(1).name = 'ksvd';
+    
+    %   Defining the parameters for KSVD
+    %   In this example we are learning 256 atoms in 20 iterations, so that
+    %   every patch in the training set can be represented with target error in
+    %   L2-norm (EData)
+    %   Type help ksvd in MATLAB prompt for more options.
+    
+    Edata=sqrt(prod(SMALL.Problem.blocksize)) * SMALL.Problem.sigma * SMALL.Problem.gain;
+    SMALL.DL(1).param=struct(...
+        'Edata', Edata,...
+        'initdict', SMALL.Problem.initdict,...
+        'dictsize', SMALL.Problem.p,...
+        'iternum', 20,...
+        'memusage', 'high');
+    
+    %   Learn the dictionary
+    
+    SMALL.DL(1) = SMALL_learn(SMALL.Problem, SMALL.DL(1));
+    
+    %   Set SMALL.Problem.A dictionary
+    %   (backward compatiblity with SPARCO: solver structure communicate
+    %   only with Problem structure, ie no direct communication between DL and
+    %   solver structures)
+    
+    SMALL.Problem.A = SMALL.DL(1).D;
+    
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    
+    
+    SMALL.solver(1)=SMALL_init_solver;
+    
+    % Defining the parameters needed for denoising
+    
+    SMALL.solver(1).toolbox='ompbox';
+    SMALL.solver(1).name='ompdenoise';
+    
+    %   Denoising the image - SMALL_denoise function is similar to SMALL_solve,
+    %   but backward compatible with KSVD definition of denoising
+    
+    SMALL.solver(1)=SMALL_denoise(SMALL.Problem, SMALL.solver(1));
+    
+    %%
+    %   Use SPAMS Online Dictionary Learning Algorithm
+    %   to Learn overcomplete dictionary (Julien Mairal 2009)
+    %   (If you have not installed SPAMS please comment the following two cells)
+    
+    %   Initialising Dictionary structure
+    %   Setting Dictionary structure fields (toolbox, name, param, D and time)
+    %   to zero values
+    
+    SMALL.DL(2)=SMALL_init_DL();
+    
+    %   Defining fields needed for dictionary learning
+    
+    SMALL.DL(2).toolbox = 'SPAMS';
+    SMALL.DL(2).name = 'mexTrainDL';
+    
+    %   Type 'help mexTrainDL in MATLAB prompt for explanation of parameters.
+    
+    SMALL.DL(2).param=struct(...
+        'D', SMALL.Problem.initdict,...
+        'K', SMALL.Problem.p,...
+        'lambda', 2,...
+        'iter', 300,...
+        'mode', 3,...
+        'modeD', 0 );
+    
+    %   Learn the dictionary
+    
+    SMALL.DL(2) = SMALL_learn(SMALL.Problem, SMALL.DL(2));
+    
+    %   Set SMALL.Problem.A dictionary
+    %   (backward compatiblity with SPARCO: solver structure communicate
+    %   only with Problem structure, ie no direct communication between DL and
+    %   solver structures)
+    
+    SMALL.Problem.A = SMALL.DL(2).D;
+    
+    
+    %%
+    %   Initialising solver structure
+    %   Setting solver structure fields (toolbox, name, param, solution,
+    %   reconstructed and time) to zero values
+    
+    SMALL.solver(2)=SMALL_init_solver;
+    
+    % Defining the parameters needed for denoising
+    
+    SMALL.solver(2).toolbox='ompbox';
+    SMALL.solver(2).name='ompdenoise';
+    
+    %   Denoising the image - SMALL_denoise function is similar to SMALL_solve,
+    %   but backward compatible with KSVD definition of denoising
+    
+    SMALL.solver(2)=SMALL_denoise(SMALL.Problem, SMALL.solver(2));
+    
+    
+    
+    %% show results %%
+    %   This will show denoised images and dictionaries for all training sets.
+    %   If you are not interested to see them and do not want clutter your
+    %   screen comment following line
+    
+    SMALL_ImgDeNoiseResult(SMALL);
+    
+    time(1,i) = SMALL.DL(1).time;
+    psnr(1,i) = SMALL.solver(1).reconstructed.psnr;
+    
+    time(2,i) = SMALL.DL(2).time;
+    psnr(2,i) = SMALL.solver(2).reconstructed.psnr;
+    
+    clear SMALL
+end
+
+%% show time and psnr %%
+figure('Name', 'KSVD vs SPAMS');
+
+subplot(1,2,1); plot(Training_size, time(1,:), 'ro-', Training_size, time(2,:), 'b*-');
+legend('KSVD','SPAMS',0);
+title('Time vs Training size');
+subplot(1,2,2); plot(Training_size, psnr(1,:), 'ro-', Training_size, psnr(2,:), 'b*-');
+legend('KSVD','SPAMS',0);
+title('PSNR vs Training size');
\ No newline at end of file
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/examples/Pierre Villars/Pierre_Villars_Example.m	Mon Mar 22 10:45:01 2010 +0000
@@ -0,0 +1,87 @@
+%%  Pierre Villars Example
+%   This example is based on the experiment suggested by Professor Pierre
+%   Vandergheynst on the SMALL meeting in Villars.
+%   The idea behind is to use patches from source image as a dictonary in
+%   which we represent target image using matching pursuit algorithm.
+%   Calling Pierre_Problem function to get src image to be used as dictionary
+%   and target image to be represented using MP with 3 paches from source image
+%
+%
+%   Ivan Damnjanovic 2010
+%%
+
+clear all;
+
+%   Defining the Problem structure
+
+SMALL.Problem = Pierre_Problem();
+
+%   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:');
+
+%   Using ten different dictionary sizes. First dictionary will contain all
+%   patches from the source image and last one will have only
+%   num_src_patches/2^9 atoms representing equidistant patches taken from
+%   the source image.
+
+n =10;
+dictsize=zeros(1,n);
+time = zeros(1,n);
+psnr = zeros(1,n);
+
+for i=1:n
+    
+    
+    %   Set reconstruction function
+    
+    SMALL.Problem.reconstruct=@(x) Pierre_reconstruct(x, SMALL.Problem);
+    
+    SMALL.solver(i)=SMALL_init_solver;
+    
+    %   Defining the parameters sparse representation
+    
+    SMALL.solver(i).toolbox='SMALL';
+    SMALL.solver(i).name='SMALL_MP';
+    
+    %   Parameters needed for matching pursuit (max number of atoms is 3
+    %   and residual error goal is 1e14
+    
+    SMALL.solver(i).param=sprintf('%d, 1e-14',3);
+    
+    % Represent the image using the source image patches as dictionary
+    
+    SMALL.solver(i)=SMALL_solve(SMALL.Problem, SMALL.solver(i));
+    
+    
+    dictsize(1,i) = size(SMALL.Problem.A,2);
+    time(1,i) = SMALL.solver(i).time;
+    psnr(1,i) = SMALL.solver(i).reconstructed.psnr;
+    
+    %   Set new SMALL.Problem.A dictionary taking every second patch from
+    %   previous dictionary
+    
+    SMALL.Problem.A=SMALL.Problem.A(:,1:2:dictsize(1,i));
+    
+    
+    %%  show reconstructed image %%
+    figure('Name', sprintf('dictsize=%d', dictsize(1,i)));
+    
+    imshow(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 ));
+    
+    
+end
+
+%%  plot time and psnr given dictionary size %%
+figure('Name', 'time and psnr');
+
+subplot(1,2,1); plot(dictsize(1,:), time(1,:), 'ro-');
+title('Time vs number of source image patches used');
+subplot(1,2,2); plot(dictsize(1,:), psnr(1,:), 'b*-');
+title('PSNR vs number of source image patches used');
\ No newline at end of file
--- a/examples/private/add_dc.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,33 +0,0 @@
-function x = add_dc(y,dc,columns)
-%ADD_DC Add DC channel to signals.
-%   X = ADD_DC(Y,DC) adds the specified DC value to the (possibly
-%   multi-dimensional) signal Y, returning the result as X. DC should be a
-%   scalar value.
-%
-%   X = ADD_DC(Y,DC,'columns') where Y is a 2D matrix and DC is an array of
-%   length size(Y,2), treats the columns of Y as individual 1D signals, 
-%   adding to each one the corresponding DC value from the DC array. X is
-%   the same size as Y and contains the resulting signals.
-%
-%   See also REMOVE_DC.
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2009
-
-
-if (nargin==3 && strcmpi(columns,'columns')), columns = 1;
-else columns = 0;
-end
-
-if (columns)
-  x = addtocols(y,dc);
-else
-  x = y + dc;
-end
-
-
-
--- a/examples/private/addtocols.c	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,85 +0,0 @@
-/**************************************************************************
- *
- * File name: addtocols.c
- *
- * Ron Rubinstein
- * Computer Science Department
- * Technion, Haifa 32000 Israel
- * ronrubin@cs
- *
- * Last Updated: 19.4.2009
- *
- *************************************************************************/
-
-
-#include "mex.h"
-
-
-/* Input Arguments */
-
-#define	X_IN	prhs[0]
-#define V_IN  prhs[1]
-
-
-/* Output Arguments */
-
-#define	Y_OUT	plhs[0]
-
-
-void mexFunction(int nlhs, mxArray *plhs[], 
-		             int nrhs, const mxArray*prhs[])
-     
-{ 
-    double *x, *y, *v, *xend;  
-    mwSize m,n,m1,n1;
-    mwIndex counter; 
-    
-    
-    /* Check for proper number of arguments */
-    
-    if (nrhs != 2) {
-      mexErrMsgTxt("Two input arguments required."); 
-    } else if (nlhs > 1) {
-      mexErrMsgTxt("Too many output arguments."); 
-    } 
-    
-    
-    /* Check the the input dimensions */ 
-    
-    m = mxGetM(X_IN);
-    n = mxGetN(X_IN);
-    if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || mxGetNumberOfDimensions(X_IN)>2) {
-      mexErrMsgTxt("ADDTOCOLS requires that X be a double matrix.");
-    }
-    m1 = mxGetM(V_IN);
-    n1 = mxGetN(V_IN);
-    if (!mxIsDouble(V_IN) || mxIsComplex(V_IN) || (m1!=1 && n1!=1)) {
-      mexErrMsgTxt("ADDTOCOLS requires that V be a double vector.");
-    }
-    if ((m1==1 && n1!=n) || (n1==1 && m1!=n)) {
-      mexErrMsgTxt("Error in ADDTOCOLS: dimensions of V and X must agree.");
-    }
-    
-    
-    /* Create a matrix for the return argument */ 
-    Y_OUT = mxCreateDoubleMatrix(m, n, mxREAL);
-         
-    
-    /* Assign pointers to the various parameters */ 
-    x = mxGetPr(X_IN); 
-    v = mxGetPr(V_IN);
-    y = mxGetPr(Y_OUT);
-            
-    
-    /* Do the actual computation */
-    
-    xend = x+(m*n);
-    counter = 0;
-    while (x<xend) {
-      (*y) = (*x) + (*v);
-      y++; x++; counter++;
-      if (counter==m) {v++; counter=0;}
-    }
-    
-    return;
-}
--- a/examples/private/addtocols.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,13 +0,0 @@
-%ADDTOCOLS Add values to the columns of a matrix.
-%  Y=ADDTOCOLS(X,V) adds to each column of the MxN matrix X the
-%  corresponding value from the N-element vector V.
-%
-%  See also NORMCOLS.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  June 2005
Binary file examples/private/addtocols.mexw32 has changed
--- a/examples/private/collincomb.c	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,165 +0,0 @@
-/**************************************************************************
- *
- * File name: collincomb.c
- *
- * Ron Rubinstein
- * Computer Science Department
- * Technion, Haifa 32000 Israel
- * ronrubin@cs
- *
- * Last Updated: 21.5.2009
- *
- *************************************************************************/
-
-
-#include "mex.h"
-
-
-/* Input Arguments */
-
-#define	A_IN	   prhs[0]
-#define ROWS_IN  prhs[1]
-#define COLS_IN1 prhs[1]
-#define COLS_IN2 prhs[2]
-#define X_IN1    prhs[2]
-#define X_IN2    prhs[3]
-
-
-/* Output Arguments */
-
-#define	Y_OUT	plhs[0]
-
-
-void mexFunction(int nlhs, mxArray *plhs[],
-int nrhs, const mxArray*prhs[])
-
-{
-  double *A, *x, *y, *rows, *cols;
-  mwSize m,n,m1,n1,m2,n2,rownum,colnum;
-  mwIndex col_id,*row_ids,i,j;
-  int rownumspecified=0;
-  
-  
-  /* Check for proper number of arguments */
-  
-  if (nrhs!=3 && nrhs!=4) {
-    mexErrMsgTxt("Invalid number of arguments.");
-  } else if (nlhs > 1) {
-    mexErrMsgTxt("Too many output arguments.");
-  }
-  
-  
-  /* Check the input dimensions */
-  
-  m = mxGetM(A_IN);
-  n = mxGetN(A_IN);
-  if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) {
-    mexErrMsgTxt("COLLINCOMB requires that A be a double matrix.");
-  }
-  
-  if (nrhs==3) {
-    
-    m1 = mxGetM(COLS_IN1);
-    n1 = mxGetN(COLS_IN1);
-    if (!mxIsDouble(COLS_IN1) || mxIsComplex(COLS_IN1) || (m1!=1 && n1!=1)) {
-      mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double.");
-    }
-    colnum = (m1 > n1) ? m1 : n1;   /* the number of columns in the linear combination */
-    
-    m2 = mxGetM(X_IN1);
-    n2 = mxGetN(X_IN1);
-    if (!mxIsDouble(X_IN1) || mxIsComplex(X_IN1) || (m2!=1 && n2!=1)) {
-      mexErrMsgTxt("COLLINCOMB requires that X be a double vector.");
-    }
-    
-    if (m2!=colnum && n2!=colnum) {
-      mexErrMsgTxt("The length of X does not match the number of columns in COLS.");
-    }
-    
-    rows = 0;
-    Y_OUT = mxCreateDoubleMatrix(m, 1, mxREAL);
-    cols = mxGetPr(COLS_IN1);
-    x = mxGetPr(X_IN1);
-  }
-  else {
-    
-    m1 = mxGetM(ROWS_IN);
-    n1 = mxGetN(ROWS_IN);
-    if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) {
-      mexErrMsgTxt("COLLINCOMB requires that ROWS be an index vector of type double.");
-    }
-    rownum = (m1 > n1) ? m1 : n1;   /* the number of rows in the linear combination */
-    rownumspecified = 1;
-    rows = mxGetPr(ROWS_IN);
-    
-    m1 = mxGetM(COLS_IN2);
-    n1 = mxGetN(COLS_IN2);
-    if (!mxIsDouble(COLS_IN2) || mxIsComplex(COLS_IN2) || (m1!=1 && n1!=1)) {
-      mexErrMsgTxt("COLLINCOMB requires that COLS be an index vector of type double.");
-    }
-    colnum = (m1 > n1) ? m1 : n1;   /* the number of columns in the linear combination */
-    
-    m2 = mxGetM(X_IN2);
-    n2 = mxGetN(X_IN2);
-    if (!mxIsDouble(X_IN2) || mxIsComplex(X_IN2) || (m2!=1 && n2!=1)) {
-      mexErrMsgTxt("COLLINCOMB requires that X be a double vector.");
-    }
-    
-    if (m2!=colnum && n2!=colnum) {
-      mexErrMsgTxt("The length of X does not match the number of columns in COLS.");
-    }
-    
-    Y_OUT = mxCreateDoubleMatrix(rownum, 1, mxREAL);
-    cols = mxGetPr(COLS_IN2);
-    x = mxGetPr(X_IN2);
-  }
-  
-  
-  /* Assign pointers to the various parameters */
-  A = mxGetPr(A_IN);
-  y = mxGetPr(Y_OUT);
-  
-  
-  if (rownumspecified) {
-    
-     /* check row indices */
-    
-    row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex));
-    
-    for (i=0; i<rownum; ++i) {
-      row_ids[i] = (mwIndex)(rows[i]+0.1)-1;
-      if (row_ids[i]<0 || row_ids[i]>=m) {
-        mexErrMsgTxt("Row index in ROWS is out of range.");
-      }
-    }
-    
-    /* Do the actual computation */
-    for (i=0; i<colnum; ++i) {
-      col_id = (mwIndex)(cols[i]+0.1)-1;
-      if (col_id<0 || col_id>=n) {
-        mexErrMsgTxt("Column index in COLS is out of range.");
-      }
-      for (j=0; j<rownum; ++j) {
-        y[j] += A[m*col_id+row_ids[j]]*x[i];
-      }
-    }
-    
-    mxFree(row_ids);
-  }
-  
-  else {
-    
-    /* Do the actual computation */
-    for (i=0; i<colnum; ++i) {
-      col_id = (mwIndex)(cols[i]+0.1)-1;
-      if (col_id<0 || col_id>=n) {
-        mexErrMsgTxt("Column index in COLS is out of range.");
-      }
-      for (j=0; j<m; ++j) {
-        y[j] += A[m*col_id+j]*x[i];
-      }
-    }
-  }
-  
-  return;
-}
--- a/examples/private/collincomb.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,27 +0,0 @@
-%COLLINCOMB Linear combination of matrix columns.
-%  Y = COLLINCOMB(A,COLS,X) computes a linear combination of the columns of
-%  the matrix A. The column indices are specified in the vector COLS, and
-%  the correspoinding coefficients are specified in the vector X. The
-%  vectors COLS and X must be of the same length. 
-%  
-%  The call Y = COLLINCOMB(A,COLS,X) is essentially equivalent to
-%
-%         Y = A(:,COLS)*X .
-%
-%  However, it is implemented much more efficiently.
-%
-%  Y = COLLINCOMB(A,ROWS,COLS,X) only works on the rows of A specified
-%  in ROWS, returning a vector of length equal to ROWS. This call is
-%  essentially equivalent to the command
-%
-%         Y = A(ROWS,COLS)*X .
-%
-%  See also ROWLINCOMB.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2009
Binary file examples/private/collincomb.mexw32 has changed
--- a/examples/private/countcover.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-function cnt = countcover(sz,blocksize,stepsize)
-%COUNTCOVER Covering of signal samples by blocks
-%  CNT = COUNTCOVER(SZ,BLOCKSIZE,STEPSIZE) assumes a p-dimensional signal
-%  of size SZ=[N1 N2 ... Np] covered by (possibly overlapping) blocks of
-%  size BLOCKSIZE=[M1 M2 ... Mp]. The blocks start at position (1,1,..,1)
-%  and are shifted between them by steps of size STEPSIZE=[S1 S2 ... Sp].
-%  COUNTCOVER returns a matrix the same size as the signal, containing in
-%  each entry the number of blocks covering that sample.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  August 2008
-
-
-cnt = ones(sz);
-for k = 1:length(sz)
-  
-  % this code is modified from function NDGRID, so it computes one
-  % output argument of NDGRID at a time (to conserve memory)
-  ids = (1:sz(k))';
-  s = sz; s(k) = [];
-  ids = reshape(ids(:,ones(1,prod(s))),[length(ids) s]);
-  ids = permute(ids,[2:k 1 k+1:length(sz)]);
-  
-  cnt = cnt .* max( min(floor((ids-1)/stepsize(k)),floor((sz(k)-blocksize(k))/stepsize(k))) - ...
-                    max(ceil((ids-blocksize(k))/stepsize(k)),0) + 1 , 0 );
-end
--- a/examples/private/diag_ids.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,30 +0,0 @@
-function id = diag_ids(x,k)
-%DIAG_IDS Indices of matrix diagonal elements.
-%  ID = DIAG_IDS(X) returns the indices of the main diagonal of X.
-%
-%  ID = DIAG_IDS(X,K) returns the indices of the K-th diagonal. K=0
-%  represents the main diagonal, positive values are above the main
-%  diagonal and negative values are below the main diagonal.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  September 2006
-
-
-if (nargin==1), k=0; end
-
-[m,n] = size(x);
-l = max(m,n);
-
-if (k<=0)
-  id = (0:l-1)*m + (1:l) - k;
-else
-  id = (0:l-1)*m + (1:l) + k*m;
-end
-
-if (l-k>m), id = id(1:end-(l-k-m)); end
-if (l+k>n), id = id(1:end-(l+k-n)); end
--- a/examples/private/dictdist.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,61 +0,0 @@
-function [dist,ratio] = dictdist(approxD,D,epsilon)
-%DICTDIST Distance between dictionaries.
-%  [DIST,RATIO] = DICTDIST(APPROXD,D) computes the distance between the
-%  approximate dictionary APPROXD and the true dictionary D, where APPROXD
-%  is NxK and D is NxM.
-%
-%  The distance between the dictionary APPROXD and a single atom A of D is
-%  defined as:
-%
-%      DIST(APPROXD,A) = min  { 1-abs(APPROXD(:,i)' * A) }
-%                         i
-%
-%  The distance between the dictionaries APPROXD and D is defined as:
-%
-%      DIST(APPROXD,D) = sum { dist(APPROXD, D(:,k)) } / M
-%                         k
-%
-%  Note that 0 <= DIST(APPROXD,D) <= 1, where 0 implies that all atoms in D
-%  appear in APPROXD, and 1 implies that the atoms of D are orthogonal to
-%  APPROXD.
-%
-%  The similarity ratio between APPROXD and D is defined as:
-%
-%      RATIO(APPROXD,D) = #(atoms in D that appear in APPROXD) / M
-%
-%  where two atoms are considered identical when DIST(A1,A2) < EPSILON with
-%  EPSILON=0.01 by default. Note that 0 <= RATIO(APPROXD,D) <= 1, where 0
-%  means APPROXD and D have no identical atoms, and 1 means that all atoms
-%  of D appear in APPROXD.
-%
-%  [DIST,RATIO] = DICTDIST(DICT1,DICT2,EPSILON) specifies a different value
-%  for EPSILON.
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  October 2007
-
-
-if (nargin < 3), epsilon = 0.01; end
-
-[n,m] = size(D);
-
-approxD = normcols(approxD*spdiag(sign(approxD(1,:))));
-D = normcols(D*spdiag(sign(D(1,:))));
-
-identical_atoms = 0;
-dist = 0;
-
-for i = 1:m
-  atom = D(:,i);
-  distances = 1-abs(atom'*approxD);
-  mindist = min(distances);
-  dist = dist + mindist;
-  identical_atoms = identical_atoms + (mindist < epsilon);
-end
-
-dist = dist / m;
-ratio = identical_atoms / m;
--- a/examples/private/imnormalize.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,19 +0,0 @@
-function y = imnormalize(x)
-%IMNORMALIZE Normalize image values.
-%  Y = IMNORMALIZE(X) linearly transforms the intensity values of the image
-%  X to tightly cover the range [0,1]. If X has more than one channel, the
-%  channels are handled as one and normalized using the same transform.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  May 2004
-
-
-maxval = max(x(:));
-minval = min(x(:));
-
-y = (x-minval) / (maxval-minval);
--- a/examples/private/iswhole.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-function z = iswhole(x,epsilon)
-%ISWHOLE Determine whole numbers with finite precision.
-%  Z = ISWHOLE(X,EPSILON) returns a matrix the same size as X, containing
-%  1's where X is whole up to precision EPSILON, and 0's elsewhere. 
-%
-%  Z = ISWHOLE(X) uses the default value of EPSILON=1e-6.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  May 2005
-
-if (nargin<2), epsilon = 1e-6; end
-z = abs(round(x)-x) < epsilon;
--- a/examples/private/make.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,40 +0,0 @@
-function make
-%MAKE Build the KSVDBox MEX support files.
-%  MAKE compiles the KSVDBox supporting MEX functions, using Matlab's
-%  default MEX compiler. If the MEX compiler has not been set-up before,
-%  please run 
-%
-%    mex -setup
-%
-%  before using this MAKE file.
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2009
-
-
-% detect platform 
-
-compstr = computer;
-is64bit = strcmp(compstr(end-1:end),'64');
-
-
-% compilation parameters
-
-compile_params = cell(0);
-if (is64bit)
-  compile_params{1} = '-largeArrayDims';
-end
-
-
-% Compile files %
-
-sourcefiles = {'addtocols.c','collincomb.c','rowlincomb.c'};
-
-for i = 1:length(sourcefiles)
-  printf('Compiling %s...', sourcefiles{i});
-  mex(sourcefiles{i},compile_params{:});
-end
--- a/examples/private/normcols.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,17 +0,0 @@
-function y = normcols(x)
-%NORMCOLS Normalize matrix columns.
-%  Y = NORMCOLS(X) normalizes the columns of X to unit length, returning
-%  the result as Y.
-%
-%  See also ADDTOCOLS.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2009
-
-
-y = x*spdiag(1./sqrt(sum(x.*x)));
--- a/examples/private/printf.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-function str = printf(varargin)
-%PRINTF Print formatted text to screen.
-%  PRINTF(FMT,VAL1,VAL2,...) formats the data in VAL1,VAL2,... according to
-%  the format string FMT, and prints the result to the screen.
-%
-%  The call to PRINTF(FMT,VAL1,VAL2,...) simply invokes the call
-%  DISP(SPRINTF(FMT,VAL1,VAL2,...)). For a complete description of the
-%  format string options see function SPRINTF.
-%
-%  STR = PRINTF(...) also returns the formatted string.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2008
-
-
-if (nargout>0)
-  str = sprintf(varargin{:});
-  disp(str);
-else
-  disp(sprintf(varargin{:}));
-end
--- a/examples/private/reggrid.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,136 +0,0 @@
-function [varargout] = reggrid(sz,num,mode)
-%REGGRID Regular sampling grid.
-%  [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM) returns the indices
-%  of a regular uniform sampling grid over a p-dimensional matrix with
-%  dimensions N1xN2x...xNp. NUM is the minimal number of required samples,
-%  and it is ensured that the actual number of samples, given by
-%  length(I1)xlength(I2)x...xlength(Ip), is at least as large as NUM.
-%
-%  [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM,'MODE') specifies the
-%  method for distributing the samples along each dimension. Valid modes
-%  include 'eqdist' (the default mode) and 'eqnum'. 'eqdist' indicates an
-%  equal distance between the samples in each dimension, while 'eqnum'
-%  indicates an equal number of samples in each dimension.
-%
-%  Notes about MODE:
-%
-%    1. The 'eqnum' mode will generally fail when the p-th root of NUM
-%    (i.e. NUM^(1/p)) is larger than min([N1 N2 ... Np]). Thus 'eqdist' is
-%    the more useful choice for sampling an arbitrary number of samples
-%    from the matrix (up to the total number of matrix entries).
-%  
-%    2. In both modes, the equality (of the distance between samples, or
-%    the number of samples in each dimension) is only approximate. This is
-%    because REGGRID attempts to maintain the appropriate equality while at
-%    the same time find a sampling pattern where the total number of
-%    samples is as close as possible to NUM. In general, the larger {Ni}
-%    and NUM are, the tighter the equality.
-%
-%  Example: Sample a set of blocks uniformly from a 2D image.
-%
-%    n = 512; blocknum = 20000; blocksize = [8 8];
-%    im = rand(n,n);
-%    [i1,i2] = reggrid(size(im)-blocksize+1, blocknum);
-%    blocks = sampgrid(im, blocksize, i1, i2);
-%
-%  See also SAMPGRID.
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  November 2007
-
-dim = length(sz);
-
-if (nargin<3)
-  mode = 'eqdist';
-end
-
-if (any(sz<1))
-  error(['Invalid matrix size : [' num2str(sz) ']']);
-end
-
-if (num > prod(sz))
-  warning(['Invalid number of samples, returning maximum number of samples.']);
-elseif (num <= 0)
-  if (num < 0)
-    warning('Invalid number of samples, assuming 0 samples.');
-  end
-  for i = 1:length(sz)
-    varargout{i} = [];
-  end
-  return;
-end
-
-
-if (strcmp(mode,'eqdist'))
-  
-  % approximate distance between samples: total volume divided by number of
-  % samples gives the average volume per sample. then, taking the p-th root
-  % gives the average distance between samples
-  d = (prod(sz)/num)^(1/dim);
-  
-  % compute the initial guess for number of samples in each dimension.
-  % then, while total number of samples is too large, decrese the number of
-  % samples by one in the dimension where the samples are the most crowded.
-  % finally, do the opposite process until just passing num, so the final
-  % number of samples is the closest to num from above.
-  
-  n = min(max(round(sz/d),1),sz);   % set n so that it saturates at 1 and sz
-  
-  active_dims = find(n>1);    % dimensions where the sample num can be reduced
-  while(prod(n)>num && ~isempty(active_dims))
-    [y,id] = min((sz(active_dims)-1)./n(active_dims));
-    n(active_dims(id)) = n(active_dims(id))-1;
-    if (n(active_dims(id)) < 2)
-      active_dims = find(n>1);
-    end
-  end
-
-  active_dims = find(n<sz);    % dimensions where the sample num can be increased
-  while(prod(n)<num && ~isempty(active_dims))
-    [y,id] = max((sz(active_dims)-1)./n(active_dims));
-    n(active_dims(id)) = n(active_dims(id))+1;
-    if (n(active_dims(id)) >= sz(active_dims(id)))
-      active_dims = find(n<sz);
-    end
-  end
-
-  for i = 1:dim
-    varargout{i} = round((1:n(i))/n(i)*sz(i));
-    varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2);
-  end
-  
-elseif (strcmp(mode,'eqnum'))
-  
-  % same idea as above
-  n = min(max( ones(size(sz)) * round(num^(1/dim)) ,1),sz);
-
-  active_dims = find(n>1);
-  while(prod(n)>num && ~isempty(active_dims))
-    [y,id] = min((sz(active_dims)-1)./n(active_dims));
-    n(active_dims(id)) = n(active_dims(id))-1;
-    if (n(active_dims(id)) < 2)
-      active_dims = find(n>1);
-    end
-  end
-  
-  active_dims = find(n<sz);
-  while(prod(n)<num && ~isempty(active_dims))
-    [y,id] = max((sz(active_dims)-1)./n(active_dims));
-    n(active_dims(id)) = n(active_dims(id))+1;
-    if (n(active_dims(id)) >= sz(active_dims(id)))
-      active_dims = find(n<sz);
-    end
-  end
-  
-  for i = 1:dim
-    varargout{i} = round((1:n(i))/n(i)*sz(i));
-    varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2);
-  end
-else
-  error('Invalid sampling mode');
-end
-
--- a/examples/private/remove_dc.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,31 +0,0 @@
-function [y,dc] = remove_dc(x,columns)
-%REMOVE_DC Remove DC channel from signals.
-%   [Y,DC] = REMOVE_DC(X) removes the DC channel (i.e. the mean) from the
-%   specified (possibly multi-dimensional) signal X. Y is the DC-free
-%   signal and is the same size as X. DC is a scalar containing the mean of
-%   the signal X.
-%
-%   [Y,DC] = REMOVE_DC(X,'columns') where X is a 2D matrix, treats the
-%   columns of X as a set of 1D signals, removing the DC channel from each
-%   one individually. Y is the same size as X and contains the DC-free
-%   signals. DC is a row vector of length size(X,2) containing the means of
-%   the signals in X.
-%
-%   See also ADD_DC.
-
-
-if (nargin==2 && strcmpi(columns,'columns')), columns = 1;
-else columns = 0;
-end
-
-if (columns)
-  dc = mean(x);
-  y = addtocols(x,-dc);
-else
-  if (ndims(x)==2)  % temporary, will remove in future
-    warning('Treating 2D matrix X as a single signal and not each column individually');
-  end
-  dc = mean(x(:));
-  y = x-dc;
-end
-
--- a/examples/private/rowlincomb.c	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,148 +0,0 @@
-/**************************************************************************
- *
- * File name: rowlincomb.c
- *
- * Ron Rubinstein
- * Computer Science Department
- * Technion, Haifa 32000 Israel
- * ronrubin@cs
- *
- * Last Updated: 21.5.2009
- *
- *************************************************************************/
-
-#include "mex.h"
-
-
-/* Input Arguments */
-
-#define	X_IN	   prhs[0]
-#define A_IN     prhs[1]
-#define ROWS_IN  prhs[2]
-#define COLS_IN  prhs[3]
-
-
-/* Output Arguments */
-
-#define	Y_OUT	plhs[0]
-
-
-void mexFunction(int nlhs, mxArray *plhs[],
-int nrhs, const mxArray*prhs[])
-
-{
-  double *A, *x, *y, *rows, *cols;
-  mwSize m,n,m1,n1,m2,n2,rownum,colnum;
-  mwIndex *row_ids,*col_ids,i,j;
-  int colnumspecified=0;
-  
-  
-  /* Check for proper number of arguments */
-  
-  if (nrhs!=3 && nrhs!=4) {
-    mexErrMsgTxt("Invalid number of input arguments.");
-  } else if (nlhs > 1) {
-    mexErrMsgTxt("Too many output arguments.");
-  }
-  
-  
-  /* Check the input dimensions */
-  
-  m = mxGetM(A_IN);
-  n = mxGetN(A_IN);
-  if (!mxIsDouble(A_IN) || mxIsComplex(A_IN) || mxGetNumberOfDimensions(A_IN)>2) {
-    mexErrMsgTxt("ROWLINCOMB requires that A be a double matrix.");
-  }
-  
-  m1 = mxGetM(ROWS_IN);
-  n1 = mxGetN(ROWS_IN);
-  if (!mxIsDouble(ROWS_IN) || mxIsComplex(ROWS_IN) || (m1!=1 && n1!=1)) {
-    mexErrMsgTxt("ROWLINCOMB requires that ROWS be an index vector of type double.");
-  }
-  rownum = (m1 > n1) ? m1 : n1;   /* the number of rows in the linear combination */
-  
-  m2 = mxGetM(X_IN);
-  n2 = mxGetN(X_IN);
-  if (!mxIsDouble(X_IN) || mxIsComplex(X_IN) || ((m2!=1) && (n2!=1))) {
-    mexErrMsgTxt("ROWLINCOMB requires that X be a double vector.");
-  }
-  
-  if (m2 != rownum && n2 != rownum) {
-    mexErrMsgTxt("The length of X does not match the number of rows in ROWS.");
-  }
-  
-  if (nrhs==4) {
-    m1 = mxGetM(COLS_IN);
-    n1 = mxGetN(COLS_IN);
-    if (!mxIsDouble(COLS_IN) || mxIsComplex(COLS_IN) || (m1!=1 && n1!=1)) {
-      mexErrMsgTxt("ROWLINCOMB requires that COLS be an index vector of type double.");
-    }
-    colnum = (m1 > n1) ? m1 : n1;   /* the number of columns */
-    colnumspecified = 1;
-    cols = mxGetPr(COLS_IN);
-    
-    Y_OUT = mxCreateDoubleMatrix(1, colnum, mxREAL);
-  }
-  else {
-    cols = 0;
-    Y_OUT = mxCreateDoubleMatrix(1, n, mxREAL);
-  }
-  
-  
-  /* Assign pointers to the various parameters */
-  A = mxGetPr(A_IN);
-  rows = mxGetPr(ROWS_IN);
-  x = mxGetPr(X_IN);
-  y = mxGetPr(Y_OUT);
-  
-  
-  /* check row indices */
-  
-  row_ids = (mwIndex*)mxMalloc(rownum*sizeof(mwIndex));
-  
-  for (i=0; i<rownum; ++i) {
-    row_ids[i] = (mwIndex)(rows[i]+0.1)-1;
-    if (row_ids[i]<0 || row_ids[i]>=m) {
-      mexErrMsgTxt("Row index in ROWS is out of range.");
-    }
-  }
-  
-  
-  
-  if (colnumspecified) {
-    
-    /* check column indices */
-    col_ids = (mwIndex*)mxMalloc(colnum*sizeof(mwIndex));
-    
-    for (i=0; i<colnum; ++i) {
-      col_ids[i] = (mwIndex)(cols[i]+0.1)-1;
-      if (col_ids[i]<0 || col_ids[i]>=n) {
-        mexErrMsgTxt("Column index in COLS is out of range.");
-      }
-    }
-    
-    /* Do the actual computation */
-    for (j=0; j<colnum; ++j) {
-      for (i=0; i<rownum; ++i) {
-        y[j] += A[m*col_ids[j]+row_ids[i]]*x[i];
-      }
-    }
-    
-    mxFree(col_ids);
-  }
-  
-  else {
-    
-    /* Do the actual computation */
-    for (j=0; j<n; ++j) {
-      for (i=0; i<rownum; ++i) {
-        y[j] += A[m*j+row_ids[i]]*x[i];
-      }
-    }
-  }
-  
-  
-  mxFree(row_ids);
-  
-  return;
-}
--- a/examples/private/rowlincomb.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,26 +0,0 @@
-%ROWLINCOMB Linear combination of matrix rows.
-%  Y = ROWLINCOMB(X,A,ROWS) computes a linear combination of the rows of
-%  the matrix A. The row indices are specified in the vector ROWS, and the
-%  correspoinding coefficients are specified in the vector X. The vectors
-%  ROWS and X must be of the same length. The call Y = ROWLINCOMB(X,A,ROWS)
-%  is essentially equivalent to the command
-%
-%         Y = X'*A(ROWS,:) .
-%
-%  However, it is implemented much more efficiently.
-%
-%  Y = ROWLINCOMB(X,A,ROWS,COLS) only works on the columns of A specified
-%  in COLS, returning a vector of length equal to COLS. This call is
-%  essentially equivalent to the command
-%
-%         Y = X'*A(ROWS,COLS) .
-%
-%  See also COLLINCOMB.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2009
Binary file examples/private/rowlincomb.mexw32 has changed
--- a/examples/private/sampgrid.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,93 +0,0 @@
-function y = sampgrid(x,blocksize,varargin)
-%SAMPGRID Sample a multi-dimensional matrix on a regular grid.
-%  Y = SAMPGRID(X,BLOCKSIZE,I1,I2,...,Ip) extracts block samples of size
-%  BLOCKSIZE from the p-dimensional matrix X, arranging the samples as the
-%  column vectors of the matrix Y. The locations of the (1,1,..,1)-th
-%  elements of each block are given in the index vectors I1,I2,..Ip. The
-%  total number of samples taken is length(I1)xlength(I2)x...xlength(Ip).
-%  BLOCKSIZE should either be a p-element vector of the form [N1,N2,...Np],
-%  or a scalar N which is shorthand for the square block size [N N ... N].
-%
-%  Example: Sample a set of blocks uniformly from a 2D image.
-%
-%    n = 512; blocknum = 20000; blocksize = [8 8];
-%    im = rand(n,n);
-%    [i1,i2] = reggrid(size(im)-blocksize+1, blocknum);
-%    blocks = sampgrid(im, blocksize, i1, i2);
-%
-%  See also REGGRID.
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  November 2007
-
-
-p = ndims(x);
-
-if (numel(blocksize)==1)
-  blocksize = ones(1,p)*blocksize;
-end
-
-n = zeros(1,p);
-for i = 1:p
-  n(i) = length(varargin{i});
-end
-
-nsamps = prod(n);
-
-% create y of the same class as x
-y = zeros(prod(blocksize),nsamps,class(x));
-
-% ids() contains the index of the current block in I1..Ip
-ids = ones(p,1);
-
-% block_ids contains the indices of the current block in X
-block_ids = cell(p,1);
-for j = 1:p
-  block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1;
-end
-
-for k = 1:nsamps
-  block = x(block_ids{:});
-  y(:,k) = block(:);
-  
-  % increment ids() and block_ids{}
-  if (k<nsamps)
-    j = 1;
-    while (ids(j) == n(j))
-      ids(j) = 1;
-      block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1;
-      j = j+1;
-    end
-    ids(j) = ids(j)+1;
-    block_ids{j} = varargin{j}(ids(j)) : varargin{j}(ids(j))+blocksize(j)-1;
-  end
-end
-
-
-% 
-% p = ndims(x);
-% 
-% n = zeros(1,p);
-% for i = 1:p
-%   n(i) = length(varargin{i});
-% end
-% 
-% nsamps = prod(n);
-% 
-% % create y of the same class as x
-% y = zeros(prod(blocksize),nsamps,class(x));
-% 
-% id = cell(p,1);
-% for k = 1:nsamps
-%   [id{:}] = ind2sub(n,k);
-%   for j = 1:p
-%     id{j} = varargin{j}(id{j}) : varargin{j}(id{j})+blocksize(j)-1;
-%   end
-%   block = x(id{:});
-%   y(:,k) = block(:);
-% end
-
--- a/examples/private/secs2hms.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,28 +0,0 @@
-function [h,m,s] = secs2hms(t)
-%SECS2HMS Convert seconds to hours, minutes and seconds.
-%  [H,M,S] = SECS2HMS(T) converts the specified number of seconds T to
-%  hours, minutes and seconds. H and M are whole numbers, and S is real.
-%
-%  Example: Estimate the remaining time of a loop
-%
-%    n = 10; tic;
-%    for i = 1:n
-%      pause(1);
-%      [h,m,s] = secs2hms( (n-i)*toc/i );
-%      printf('estimated remaining time: %02d:%02d:%05.2f',h,m,s);
-%    end
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  April 2008
-
-
-s = t;
-h = fix(s/3600);
-s = rem(s,3600);
-m = fix(s/60);
-s = rem(s,60);
--- a/examples/private/spdiag.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,38 +0,0 @@
-function Y = spdiag(V,K)
-%SPDIAG Sparse diagonal matrices.
-%   SPDIAG(V,K) when V is a vector with N components is a sparse square
-%   matrix of order N+ABS(K) with the elements of V on the K-th diagonal. 
-%   K = 0 is the main diagonal, K > 0 is above the main diagonal and K < 0
-%   is below the main diagonal. 
-%
-%   SPDIAG(V) is the same as SPDIAG(V,0) and puts V on the main diagonal.
-%
-%   See also DIAG, SPDIAGS.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  June 2008
-
-
-if (nargin<2)
-  K = 0;
-end
-
-n = length(V) + abs(K);
-
-if (K>0)
-  i = 1:length(V);
-  j = K+1:n;
-elseif (K<0)
-  i = -K+1:n;
-  j = 1:length(V);
-else
-  i = 1:n;
-  j = 1:n;
-end
-
-Y = sparse(i,j,V(:),n,n);
\ No newline at end of file
--- a/examples/private/timerclear.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,37 +0,0 @@
-function timerclear()
-%TIMERCLEAR Clear all timers.
-%   TIMERCLEAR clears all currenly registered timers, invalidating all
-%   timer ids.
-%
-%   Note: since registered timers do not consume CPU power except for when
-%   the TIMER<*> functions are called, this function is only useful in
-%   situations where a large number of timers have been initialized, and
-%   there is a need to reclaim memory.
-%
-%   See also TIMERINIT, TIMERETA.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  June 2008
-
-
-global utiltbx_timer_start_times    % start times
-global utiltbx_time_lastdisp        % last display times
-global utiltbx_timer_iternums       % iteration numbers
-global utiltbx_timer_lastiter       % last queried iteration numbers
-global utiltbx_timer_name           % timer names
-global utiltbx_timer_callfun        % timer calling functions
-
-
-% clear all timers %
-
-utiltbx_timer_start_times = [];
-utiltbx_time_lastdisp = [];
-utiltbx_timer_iternums = [];
-utiltbx_timer_lastiter = [];
-utiltbx_timer_name = [];
-utiltbx_timer_callfun = [];
--- a/examples/private/timereta.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,98 +0,0 @@
-function varargout = timereta(tid,iter,delay)
-%TIMERETA Estimated remaining time.
-%   S = TIMERETA(TID,ITER) returns the estimated remaining time (in
-%   seconds) for the process associated with timer TID, assuming the
-%   process has completed ITER iterations. Note: the function will exit
-%   with an error if the timer TID does not have an associated number of
-%   iterations (see function TIMERINIT).
-%
-%   [H,M,S] = TIMERETA(TID,ITER) returns the estimated remaining time in
-%   hours, minutes and seconds.
-%
-%   TIMERETA(TID,ITER), with no output arguments, prints the estimated
-%   remaining time to the screen. The time is displayed in the format
-%
-%     TIMERNAME: iteration ITER / ITERNUM, estimated remaining time: HH:MM:SS.SS
-%
-%   If the timer has no assigned name, the display format changes to
-%
-%     Iteration ITER / ITERNUM, estimated remaining time: HH:MM:SS.SS
-%
-%   TIMERETA(TID,ITER,DELAY) only displays the remaining time if the
-%   time elapsed since the previous printout is at least DELAY seconds.
-%
-%   See also TIMERINIT, TIMERCLEAR.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  June 2008
-
-
-global utiltbx_timer_start_times 
-global utiltbx_timer_iternums
-global utiltbx_timer_lastiter
-global utiltbx_time_lastdisp
-global utiltbx_timer_name
-
-
-if (tid<1 || tid>length(utiltbx_timer_iternums))
-  error('Unknown timer id');
-end
-
-if (utiltbx_timer_iternums(tid) < 0)
-  error('Specified timer does not have an associated number of iterations');
-end
-
-% update last reported iteration number
-utiltbx_timer_lastiter(tid) = iter;
-
-% compute elapsed time
-starttime = utiltbx_timer_start_times(tid,:);
-currtime = clock;
-timediff = etime(currtime, starttime);
-
-% total iteration number
-iternum = utiltbx_timer_iternums(tid);
-
-% compute eta
-timeremain = (iternum-iter)*timediff/iter;
-
-% return eta in seconds
-if (nargout==1)
-  varargout{1} = timeremain;
-  
-% return eta in hms
-elseif (nargout==3)
-  [varargout{1}, varargout{2}, varargout{3}] = secs2hms(timeremain);
-  
-  
-% print eta
-elseif (nargout==0)
-  
-  % check last display time
-  lastdisptime = utiltbx_time_lastdisp(tid,:);
-  if (nargin>2 && etime(currtime,lastdisptime) < delay)
-    return;
-  end
-  
-  % update last display time
-  utiltbx_time_lastdisp(tid,:) = currtime;
-
-  % display timer
-  [hrs,mins,secs] = secs2hms(timeremain);
-  if (isempty(utiltbx_timer_name{tid}))
-    printf('Iteration %d / %d, estimated remaining time: %02d:%02d:%05.2f', iter, iternum, hrs, mins, secs);
-  else
-    timername = utiltbx_timer_name{tid};
-    printf('%s: iteration %d / %d, estimated remaining time: %02d:%02d:%05.2f', timername, iter, iternum, hrs, mins, secs);
-  end
-   
-% invalid number of outputs
-else
-  error('Invalid number of output arguments');
-end
-
--- a/examples/private/timerinit.m	Mon Mar 22 10:43:01 2010 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,110 +0,0 @@
-function tid = timerinit(par1,par2)
-%TIMERINIT Initialize a new timer.
-%   TID = TIMERINIT() initializes a new timer for counting elapsed time,
-%   and returns its id.
-%
-%   TID = TIMERINIT('TIMERNAME') sets the timer name to the specified
-%   string for display purposes.
-%
-%   TID = TIMERINIT(ITERNUM) initializes a new ETA timer for a process with
-%   ITERNUM iterations. An ETA timer can be used for both counting elapsed
-%   time and estimating remaining time.
-%
-%   TID = TIMERINIT('TIMERNAME',ITERNUM) sets the ETA timer name to the
-%   specified string for display purposes.
-%
-%   Example:
-%
-%     tid = timerinit(100); 
-%     for i = 1:100
-%       pause(0.07);
-%       timereta(tid,i,1);
-%     end
-%     timereta(tid,i);
-%
-%   See also TIMERETA, TIMERCLEAR.
-
-
-%  Ron Rubinstein
-%  Computer Science Department
-%  Technion, Haifa 32000 Israel
-%  ronrubin@cs
-%
-%  June 2008
-
-
-global utiltbx_timer_start_times    % start times
-global utiltbx_time_lastdisp        % last display times
-global utiltbx_timer_iternums       % iteration numbers
-global utiltbx_timer_lastiter       % last queried iteration numbers
-global utiltbx_timer_name           % timer names
-global utiltbx_timer_callfun        % timer calling functions
-
-
-% parse function arguments %
-
-if (nargin==0)
-
-  iternum = -1;
-  timername = '';
-
-elseif (nargin==1)
-
-  if (ischar(par1))
-    iternum = -1;
-    timername = par1;
-
-  elseif (isnumeric(par1) && numel(par1)==1 && par1>0)
-    iternum = par1;
-    timername = '';
-
-  else
-    error('Invalid number of iterations');
-  end
-
-elseif (nargin==2)
-
-  if (ischar(par1) && isnumeric(par2))
-    if (numel(par2)==1 && par2>0)
-      timername = par1;
-      iternum = par2;
-    else
-      error('Invalid number of iterations');
-    end
-  else
-    error('Invalid function syntax');
-  end
-
-else
-  error('Too many function parameters');
-end
-
-
-% register the timer %
-
-if (isempty(utiltbx_timer_start_times))
-  utiltbx_timer_start_times = clock;
-  utiltbx_time_lastdisp = utiltbx_timer_start_times;
-  utiltbx_timer_iternums = double(iternum);
-  utiltbx_timer_lastiter = 0;
-  utiltbx_timer_name = { timername };
-  utiltbx_timer_callfun = {};
-  tid = 1;
-else
-  utiltbx_timer_start_times(end+1,:) = clock;
-  utiltbx_time_lastdisp(end+1,:) = utiltbx_timer_start_times(end,:);
-  utiltbx_timer_iternums(end+1) = double(iternum);
-  utiltbx_timer_lastiter(end+1) = 0;
-  utiltbx_timer_name{end+1} = timername;
-  tid = size(utiltbx_timer_start_times,1);
-end
-
-
-% detect timer calling function %
-
-st = dbstack;
-if (length(dbstack) >= 2)
-  utiltbx_timer_callfun{end+1} = st(2).name;
-else
-  utiltbx_timer_callfun{end+1} = '';
-end