Mercurial > hg > smallbox
changeset 10:207a6ae9a76f version1.0
(none)
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/Pierre_Problem.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,112 @@ +function data=Pierre_Problem(src, trg, blocksize, dictsize); +%%% Generate Pierre Problem +% Ivan Damnjanovic 2010 +% +% Pierre_Problem is a part of the SMALLbox and generates the problem +% suggested by Professor Pierre Vandergheynst on the SMALL meeting in +% Villars. +% The function takes as an input: +% - src - source image matrix (if not present function promts user for +% an image file) , +% - trg - target image matrix (if not present function promts user for +% an image file) , +% - blocksize - block (patch) vertical/horizontal dimension (default 8), +% - dictsize - dictionary size (default - all patches from target +% image). +% +% The output of the function is stucture with following fields: +% - srcname - source image name, +% - imageSrc - source image matrix, +% - trgname - target image name, +% - imageTrg - Target image matrix, +% - A - dictonary with patches from the source image, +% - b - measurement matrix (i.e. patches from target image to be +% represented in dictionary A, +% - m - size of patches (default 25), +% - n - number of patches to be represented, +% - p - dictionary size, +% - blocksize - block size (default [5 5]), +% - maxval - maximum value (default - 255) +% - sparse - if 1 SMALL_solve will keep solution matrix in sparse form, +% due to memory constrains. + +%% prompt user for images %% + +%ask for source file name + +TMPpath=pwd; +FS=filesep; +if ~ exist( 'src', 'var' ) || isempty(src) +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'images']); +[filename,pathname] = uigetfile({'*.png;'},'Select a source image'); +[pathstr, name, ext, versn] = fileparts(filename); +data.srcname=name; +src = imread(filename); +src = double(src); +end; + +%ask for target file name + +if ~ exist( 'trg', 'var' ) || isempty(trg) +[filename,pathname] = uigetfile({'*.png;'},'Select a target image'); +[pathstr, name, ext, versn] = fileparts(filename); +data.trgname=name; +trg = imread(filename); +trg = double(trg); +end; +cd(TMPpath); + +%% set parameters %% + +maxval = 255; +if ~ exist( 'blocksize', 'var' ) || isempty(blocksize),blocksize = 5;end + +if ~ exist( 'dictsize', 'var' ) || isempty(dictsize), + dictsize = (size(src,1)-blocksize+1)*(size(src,2)-blocksize+1); + patch_idx=1:dictsize; +else + num_blocks_src=(size(src,1)-blocksize+1)*(size(src,2)-blocksize+1); + patch_idx=1:floor(num_blocks_src/dictsize):dictsize*floor(num_blocks_src/dictsize); +end + +p = ndims(src); +if (p==2 && any(size(src)==1) && length(blocksize)==1) + p = 1; +end + + +% blocksize % +if (numel(blocksize)==1) + blocksize = ones(1,p)*blocksize; +end +%% +%% create dictionary data %% + +S=im2col(src,blocksize,'sliding'); + +for j= 1:size(S,2) + S(:,j)=S(:,j)./norm(S(:,j)); +end + +%% create measurement matrix %% + +T=im2col(trg,blocksize, 'distinct'); + +%% output structure %% + +data.imageSrc = src; +data.imageTrg = trg; +data.A = S(:,patch_idx); +data.b = T; +data.m = size(T,1); +data.n = size(T,2); +data.p = size(data.A,2); +data.blocksize=blocksize; +data.maxval=maxval; + +% geting around out of memory problem when converting big matrix from +% sparse to full... (check SMALL_solve function) +data.sparse=1; + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/generateAMT_Learning_Problem.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,84 @@ +function data = generateAMT_Learning_Problem(nfft, windowSize, overlap) +%%% Generate Automatic Music Transcription Problem +% Ivan Damnjanovic 2010 +% +% +% generateAMT_Learning_Problem is a part of the SMALLbox and generates +% a problem that can be used for comparison of Dictionary Learning/Sparse +% Representation techniques in automatic music transcription scenario. +% The function prompts a user for an audio file (mid, wav, mat) reads it +% and generates a spectrogram given fft size (default nfft=4096), analysis +% window size (windowSize=2822), and analysis window overlap (overlap = +% 0.5). +% +% The output of the function is stucture with following fields: +% b - matrix with magnitudes of the spectrogram +% f - vector of frequencies at wihch spectrogram is computed +% windowSize - analysis window size +% overlap - analysis window overlap +% fs - sampling frequency +% m - number of frequenciy points in spectrogram +% n - number of time points in the spectrogram +% p - number of dictionary elements to be learned (eg 88 for piano) +% notesOriginal - notes of the original audio to be used for +% comparison (if midi of the original exists) +% name - name of the audio file to transcribe + +%% +FS=filesep; +if ~ exist( 'nfft', 'var' ) || isempty(nfft), nfft = 4096; end +if ~ exist( 'windowSize', 'var' ) || isempty(windowSize), windowSize = 2822; end +if ~ exist( 'overlap', 'var' ) || isempty(overlap), overlap = 0.5; end + +%% +%ask for file name +TMPpath=pwd; +[pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); +cd([pathstr1,FS,'data',FS,'audio']); +[filename,pathname] = uigetfile({'*.mat; *.mid; *.wav'},'Select a file to transcribe'); +[pathstr, name, ext, versn] = fileparts(filename); +data.name=name; + +data.notesOriginal=[]; + +if strcmp(ext,'.mid') + midi=readmidi(filename); + data.notesOriginal=midiInfo(midi); + y=midi2audio(midi); + wavwrite(y, 44100, 16, 'temp.wav'); + [x.signal, x.fs, x.nbits]=wavread('temp.wav'); + delete('temp.wav'); +elseif strcmp(ext,'.wav') + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']); + filename1=[name, '.mid']; + if exist(filename1, 'file') + midi=readmidi(filename1); + data.notesOriginal=midiInfo(midi); + end + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'wav']); + [x.signal, x.fs, x.nbits]=wavread(filename); +else + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'midi']); + filename1=[name, '.mid']; + if exist(filename1, 'file') + midi=readmidi(filename1); + data.notesOriginal=midiInfo(midi); + end + cd([pathstr1,FS, 'data', FS, 'audio', FS, 'mat']); + x=load([pathname,filename]); +end +%% +[X, frX]=spectrogram(x.signal, hanning(windowSize), overlap*windowSize, nfft, x.fs); +%% +data.b=abs(X); +data.f=frX; +data.windowSize=windowSize; +data.overlap=overlap; +data.fs=x.fs; +data.m=size(X,1); +data.n=size(X,2); + +data.p=88; %number of dictionary elements (ie notes to recover) +cd(TMPpath); + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/generateAudioDenoiseProblem.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,159 @@ +function data=generateAudioDenoiseProblem(au, trainnum, blocksize, dictsize, overlap, sigma, gain, maxval, initdict); + +% Ivan Damnjanovic 2010 +% +% +% generateAudioDenoiseProblem is part of the SMALLbox and generate a +% problem for comaprison of Dictionary Learning/Sparse Representation +% techniques in audio denoising scenario. It is based on KSVD image +% denoise demo by Ron Rubinstein (see bellow). +% The fuction takes as an optional input +% au - audio samples to be denoised +% trainnum - number of frames for training +% blocksize - 1D frame size (eg 512) +% dictsize - number of atoms to be trained +% overlap - ammount of overlaping frames between 0 and 1 +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% +% August 2009 + + +disp(' '); +disp(' ********** Denoising Problem **********'); +disp(' '); +disp(' This function reads an audio, adds random Gaussian noise,'); +disp(' that can be later denoised by using dictionary learning techniques.'); +disp(' '); + +FS=filesep; +if ~ exist( 'sigma', 'var' ) || isempty(sigma), sigma = 26.74; end +if ~ exist( 'gain', 'var' ) || isempty(gain), gain = 1.15; end + +if ~ exist( 'initdict', 'var' ) || isempty(initdict), initdict = 'odct'; end +if ~ exist( 'overlap', 'var' ) || isempty(overlap), overlap = 15/16; end +%% prompt user for wav file %% +%ask for file name + +TMPpath=pwd; +if ~ exist( 'au', 'var' ) || isempty(au) + [pathstr1, name, ext, versn] = fileparts(which('SMALLboxSetup.m')); + cd([pathstr1,FS,'data',FS,'audio',FS,'wav']); + [filename,pathname] = uigetfile({'*.wav;'},'Select a wav file'); + [pathstr, name, ext, versn] = fileparts(filename); + data.name=name; + + au = wavread(filename); + au = mean(au,2); % turn it into mono. +end; +if ~ exist( 'maxval', 'var' ) || isempty(maxval), maxval = max(au); end + +%% generate noisy audio %% + +disp(' '); +disp('Generating noisy audio...'); +sigma = max(au)/10^(sigma/20); +n = randn(size(au)) .* sigma; +aunoise = au + n;% here we can load noise audio if available + % for example: wavread('icassp06_x.wav');% + + + +%% set parameters %% + +x = aunoise; +if ~ exist( 'blocksize', 'var' ) || isempty(blocksize),blocksize = 512;end +if ~ exist( 'dictsize', 'var' ) || isempty(dictsize), dictsize = 2048;end + +if ~ exist( 'trainnum', 'var' ) || isempty(trainnum),trainnum = (size(x,1)-blocksize+1);end + + + + + +p=1; + + +% +% msgdelta = 5; +% +% verbose = 't'; +% if (msgdelta <= 0) +% verbose=''; +% msgdelta = -1; +% end +% +% +% % initial dictionary % +% +if (strcmpi(initdict,'odct')) + initdict = odctndict(blocksize,dictsize,p); +elseif (strcmpi(initdict,'data')) + clear initdict; % causes initialization using random examples +else + error('Invalid initial dictionary specified.'); +end + +if exist( 'initdict', 'var' ) + initdict = initdict(:,1:dictsize); +end + + +% noise mode % +% if (isfield(params,'noisemode')) +% switch lower(params.noisemode) +% case 'psnr' +% sigma = maxval / 10^(params.psnr/20); +% case 'sigma' +% sigma = params.sigma; +% otherwise +% error('Invalid noise mode specified'); +% end +% elseif (isfield(params,'sigma')) +% sigma = params.sigma; +% elseif (isfield(params,'psnr')) +% sigma = maxval / 10^(params.psnr/20); +% else +% error('Noise strength not specified'); +% end + +% params.Edata = sqrt(prod(blocksize)) * sigma * gain; % target error for omp +% params.codemode = 'error'; +% +% params.sigma = sigma; +% params.noisemode = 'sigma'; +% +% +% % make sure test data is not present in params +% if (isfield(params,'testdata')) +% params = rmfield(params,'testdata'); +% end + + +%%%% create training data %%% + + +X = buffer( x(1:trainnum),blocksize, overlap*blocksize); + +% remove dc in blocks to conserve memory % +% bsize = 2000; +% for i = 1:bsize:size(X,2) +% blockids = i : min(i+bsize-1,size(X,2)); +% X(:,blockids) = remove_dc(X(:,blockids),'columns'); +% end +data.Original = au; +data.Noisy = aunoise; +data.b = X; +data.m = size(X,1); +data.n = size(X,2); +data.p = dictsize; +data.blocksize=blocksize; +data.sigma = sigma; +data.gain = gain; +data.maxval = maxval; +data.initdict= initdict; +data.signalDim=1; +cd(TMPpath); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/generateImageDenoiseProblem.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,144 @@ +function data=generateImageDenoiseProblem(im, trainnum, blocksize, dictsize, sigma, gain, maxval, initdict); +%%% Generate Image Denoising Problem +% Ivan Damnjanovic 2010 +% +% generateImageDenoiseProblem is a part of the SMALLbox and generates +% a problem that can be used for comparison of Dictionary Learning/Sparse +% Representation techniques in image denoising scenario. +% The function takes as an input: +% - im - image matrix (if not present function promts user for an +% image file) , +% - trainnum - number of training samples (default - 40000) +% - blocksize - block (patch) vertical/horizontal dimension (default 8), +% - dictsize - dictionary size (default - 256), +% - sigma - noise level (default - 20), +% - noise gain (default - 1.15), +% - maxval - maximum value (default - 255) +% - initdict - initial dictionary (default - 4x overcomlete dct) +% +% The output of the function is stucture with following fields: +% - name - name of the original image (if image is read inside of the +% function) +% - Original - original image matrix, +% - Noisy - image with added noise, +% - b - training patches, +% - m - size of training patches (default 64), +% - n - number of training patches, +% - p - number of dictionary elements to be learned, +% - blocksize - block size (default [8 8]), +% - sigma - noise level, +% - noise gain (default - 1.15), +% - maxval - maximum value (default - 255) +% - initdict - initial dictionary (default - 4x overcomlete dct) +% - signalDim - signal dimension (default - 2) +% +% Based on KSVD denoise demo by Ron Rubinstein +% See also KSVDDENOISEDEMO and KSVDDEMO. +% Ron Rubinstein +% Computer Science Department +% Technion, Haifa 32000 Israel +% ronrubin@cs +% August 2009 +%% +disp(' '); +disp(' ********** Denoising Problem **********'); +disp(' '); +disp(' This function reads an image, adds random Gaussian noise,'); +disp(' that can be later denoised by using dictionary learning techniques.'); +disp(' '); + + +%% prompt user for image %% +%ask for file name +FS=filesep; +TMPpath=pwd; +if ~ exist( 'im', 'var' ) || isempty(im) + [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); + data.name=name; + im = imread(filename); + im = double(im); +end; +cd(TMPpath); + +%% check input parameters %% + +if ~ exist( 'blocksize', 'var' ) || isempty(blocksize),blocksize = 8;end +if ~ exist( 'dictsize', 'var' ) || isempty(dictsize), dictsize = 256;end +if ~ exist( 'trainnum', 'var' ) || isempty(trainnum),trainnum = 40000;end +if ~ exist( 'sigma', 'var' ) || isempty(sigma), sigma = 20; end +if ~ exist( 'gain', 'var' ) || isempty(gain), gain = 1.15; end +if ~ exist( 'maxval', 'var' ) || isempty(maxval), maxval = 255; end +if ~ exist( 'initdict', 'var' ) || isempty(initdict), initdict = 'odct'; end + +%% generate noisy image %% + +disp(' '); +disp('Generating noisy image...'); + +n = randn(size(im)) * sigma; +imnoise = im + n; + +%% set parameters %% + +x = imnoise; +p = ndims(x); + +if (p==2 && any(size(x)==1) && length(blocksize)==1) + p = 1; +end + +% blocksize % + +if (numel(blocksize)==1) + blocksize = ones(1,p)*blocksize; +end + +if (strcmpi(initdict,'odct')) + initdict = odctndict(blocksize,dictsize,p); +elseif (strcmpi(initdict,'data')) + clear initdict; % causes initialization using random examples +else + error('Invalid initial dictionary specified.'); +end + +if exist( 'initdict', 'var' ) + initdict = initdict(:,1:dictsize); +end + +%%%% create training data %%% + +ids = cell(p,1); +if (p==1) + ids{1} = reggrid(length(x)-blocksize+1, trainnum, 'eqdist'); +else + [ids{:}] = reggrid(size(x)-blocksize+1, trainnum, 'eqdist'); +end +X = sampgrid(x,blocksize,ids{:}); + +% remove dc in blocks to conserve memory % + +bsize = 2000; +for i = 1:bsize:size(X,2) + blockids = i : min(i+bsize-1,size(X,2)); + X(:,blockids) = remove_dc(X(:,blockids),'columns'); +end + +%% output structure %% + +data.Original = im; +data.Noisy = imnoise; +data.b = X; +data.m = size(X,1); +data.n = size(X,2); +data.p = dictsize; +data.blocksize=blocksize; +data.sigma = sigma; +data.gain = gain; +data.maxval = maxval; +data.initdict= initdict; +data.signalDim=2; + +end %% end of function \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/add_dc.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,33 @@ +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 + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/addtocols.c Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,85 @@ +/************************************************************************** + * + * 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; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/addtocols.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,13 @@ +%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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/collincomb.c Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,165 @@ +/************************************************************************** + * + * 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; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/collincomb.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,27 @@ +%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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/countcover.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,31 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/diag_ids.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,30 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/dictdist.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,61 @@ +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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/imnormalize.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,19 @@ +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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/iswhole.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,17 @@ +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;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/make.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,40 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/normcols.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,17 @@ +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)));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/printf.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,26 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/reggrid.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,136 @@ +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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/remove_dc.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,31 @@ +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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/rowlincomb.c Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,148 @@ +/************************************************************************** + * + * 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; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/rowlincomb.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,26 @@ +%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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/sampgrid.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,93 @@ +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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/secs2hms.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,28 @@ +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);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/spdiag.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,38 @@ +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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/timerclear.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,37 @@ +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 = [];
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/timereta.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,98 @@ +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 +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Problems/private/timerinit.m Mon Mar 22 15:06:25 2010 +0000 @@ -0,0 +1,110 @@ +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