Mercurial > hg > smallbox
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
--- 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
--- 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
--- 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