Mercurial > hg > smallbox
changeset 65:55faa9b5d1ac
(none)
author | idamnjanovic |
---|---|
date | Wed, 16 Mar 2011 13:41:02 +0000 |
parents | 8288a23f041f |
children | a02503d91c8d |
files | DL/RLS-DLA/SMALL_rlsdla 05032011.m DL/RLS-DLA/SMALL_rlsdla.m DL/RLS-DLA/SMALL_rlsdla1.m DL/RLS-DLA/SMALL_rlsdlaFirstClustTry.m Problems/generateImageDenoiseProblem.m To Do list.m examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m |
diffstat | 7 files changed, 13 insertions(+), 1466 deletions(-) [+] |
line wrap: on
line diff
--- a/DL/RLS-DLA/SMALL_rlsdla 05032011.m Tue Mar 15 15:30:54 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,488 +0,0 @@ -function Dictionary = SMALL_rlsdla(X, params) - - - - - -global CODE_SPARSITY CODE_ERROR codemode -global MEM_LOW MEM_NORMAL MEM_HIGH memusage -global ompfunc ompparams exactsvd - -CODE_SPARSITY = 1; -CODE_ERROR = 2; - -MEM_LOW = 1; -MEM_NORMAL = 2; -MEM_HIGH = 3; - - -% p = randperm(size(X,2)); - - % coding mode % -%X_norm=sqrt(sum(X.^2, 1)); -% X_norm_1=sum(abs(X)); -% X_norm_inf=max(abs(X)); -%[X_norm_sort, p]=sort(X_norm);%, 'descend'); -% [X_norm_sort1, p5]=sort(X_norm_1);%, 'descend'); - -% if (isfield(params,'codemode')) -% switch lower(params.codemode) -% case 'sparsity' -% codemode = CODE_SPARSITY; -% thresh = params.Tdata; -% case 'error' -% codemode = CODE_ERROR; -% thresh = params.Edata; -% otherwise -% error('Invalid coding mode specified'); -% end -% elseif (isfield(params,'Tdata')) -% codemode = CODE_SPARSITY; -% thresh = params.Tdata; -% elseif (isfield(params,'Edata')) -% codemode = CODE_ERROR; -% thresh = params.Edata; -% -% else -% error('Data sparse-coding target not specified'); -% end - -thresh = params.Edata; -% max number of atoms % - -% if (codemode==CODE_ERROR && isfield(params,'maxatoms')) -% ompparams{end+1} = 'maxatoms'; -% ompparams{end+1} = params.maxatoms; -% end - - -% memory usage % - -if (isfield(params,'memusage')) - switch lower(params.memusage) - case 'low' - memusage = MEM_LOW; - case 'normal' - memusage = MEM_NORMAL; - case 'high' - memusage = MEM_HIGH; - otherwise - error('Invalid memory usage mode'); - end -else - memusage = MEM_NORMAL; -end - - -% iteration count % - -if (isfield(params,'iternum')) - iternum = params.iternum; -else - iternum = 10; -end - - -% omp function % - -if (codemode == CODE_SPARSITY) - ompfunc = @omp; -else - ompfunc = @omp2; -end - - -% % status messages % -% -% printiter = 0; -% printreplaced = 0; -% printerr = 0; -% printgerr = 0; -% -% verbose = 't'; -% msgdelta = -1; -% - -% -% for i = 1:length(verbose) -% switch lower(verbose(i)) -% case 'i' -% printiter = 1; -% case 'r' -% printiter = 1; -% printreplaced = 1; -% case 't' -% printiter = 1; -% printerr = 1; -% if (isfield(params,'testdata')) -% printgerr = 1; -% end -% end -% end -% -% if (msgdelta<=0 || isempty(verbose)) -% msgdelta = -1; -% end -% -% ompparams{end+1} = 'messages'; -% ompparams{end+1} = msgdelta; -% -% -% -% % compute error flag % -% -% comperr = (nargout>=3 || printerr); -% -% -% % validation flag % -% -% testgen = 0; -% if (isfield(params,'testdata')) -% testdata = params.testdata; -% if (nargout>=4 || printgerr) -% testgen = 1; -% end -% end - -% -% % data norms % -% -% XtX = []; XtXg = []; -% if (codemode==CODE_ERROR && memusage==MEM_HIGH) -% XtX = colnorms_squared(data); -% if (testgen) -% XtXg = colnorms_squared(testdata); -% end -% end - - -% mutual incoherence limit % - -if (isfield(params,'muthresh')) - muthresh = params.muthresh; -else - muthresh = 0.99; -end -if (muthresh < 0) - error('invalid muthresh value, must be non-negative'); -end - - - - - -% determine dictionary size % - -if (isfield(params,'initdict')) - if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:)))) - dictsize = length(params.initdict); - else - dictsize = size(params.initdict,2); - end -end -if (isfield(params,'dictsize')) % this superceedes the size determined by initdict - dictsize = params.dictsize; -end - -if (size(X,2) < dictsize) - error('Number of training signals is smaller than number of atoms to train'); -end - - -% initialize the dictionary % - -if (isfield(params,'initdict')) - if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:)))) - D = X(:,params.initdict(1:dictsize)); - else - if (size(params.initdict,1)~=size(X,1) || size(params.initdict,2)<dictsize) - error('Invalid initial dictionary'); - end - D = params.initdict(:,1:dictsize); - end -else - data_ids = find(colnorms_squared(X) > 1e-6); % ensure no zero data elements are chosen - perm = randperm(length(data_ids)); - D = X(:,data_ids(perm(1:dictsize))); -end - - -% normalize the dictionary % - -% D = normcols(D); -% DtD=D'*D; - -err = zeros(1,iternum); -gerr = zeros(1,iternum); - -if (codemode == CODE_SPARSITY) - errstr = 'RMSE'; -else - errstr = 'mean atomnum'; -end -%X(:,p(X_norm_sort<thresh))=0; -% if (iternum==4) -% X_im=col2imstep(X, [256 256], [8 8]); -% else -% X_im=col2imstep(X, [512 512], [8 8]); -% end -% figure(10); imshow(X_im); - -%p1=p(cumsum(X_norm_sort)./[1:size(X_norm_sort,2)]>thresh); -%p1=p(X_norm_sort>thresh); -% X(:,setxor(p1,1:end))=0; -% X_im=col2imstep(X, [256 256], [8 8]); -% figure(10); imshow(X_im); -% if iternum==2 -% D(:,1)=D(:,2); -% end -%p1=p1(p2(1:40000)); -%end-min(40000, end)+1:end));%1:min(40000, end))); -%p1 = randperm(size(data,2));%size(data,2) -%data=data(:,p1); - -C=(100000*thresh)*eye(dictsize); -% figure(11); -w=zeros(dictsize,1); -replaced=zeros(dictsize,1); -u=zeros(dictsize,1); -% dictimg = showdict(D,[8 8],round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); -% figure(11);imshow(imresize(dictimg,2,'nearest')); -% pause(1); -lambda=0.9997%0.99986;%3+0.0001*params.linc; -for j=1:1 - %data=X; -if size(X,2)>40000 - p2 = randperm(size(X,2)); - - p2=sort(p2(1:40000));%min(floor(size(p1,2)/2),40000))); - size(p2,2) - data=X(:,p2); -elseif size(X,2)>0 - %p2 = randperm(size(p1,2)); - size(X,2) - data=X; -else - break; -end -% figure(1); -% plot(sqrt(sum(data.^2, 1))); -% a=size(data,2)/4; -% lambda0=0.99;%1-16/numS+iternum*0.0001-0.0002 -%C(1,1)=0; -modi=1000; -for i = 1:size(data,2) -% if norm(data(:,i))>thresh - % par.multA= @(x,par) multMatr(D,x); % user function y=Ax - % par.multAt=@(x,par) multMatrAdj(D,x); % user function y=A'*x - % par.y=data(:,i); - % w=SolveFISTA(D,data(:,i),'lambda',0.5*thresh); - % w=sesoptn(zeros(dictsize,1),par.func_u, par.func_x, par.multA, par.multAt,options,par); - %w = SMALL_chol(D,data(:,i), 256,32, thresh);% - %w = sparsecode(data(:,i), D, [], [], thresh); - w = omp2mex(D,data(:,i),[],[],[],thresh,0,-1,-1,0); - - %w(find(w<1))=0; - %^2; -% lambda(i)=1-0.001/(1+i/a); -% if i<a -% lambda(i)=1-0.001*(1-(i/a)); -% else -% lambda(i)=1; -% end -% param.lambda=thresh; -% param.mode=2; -% param.L=32; -% w=mexLasso(data(:,i), D, param); - spind=find(w); - %replaced(spind)=replaced(spind)+1; - %-0.001*(1/2)^(i/a); -% w_sp(i)=nnz(w); - residual = data(:,i) - D * w; - %if ~isempty(spind) - %i - if (j==1) - C = C *(1/ lambda); - end - u = C(:,spind) * w(spind); - - %spindu=find(u); - % v = D' * residual; - - alfa = 1/(1 + w' * u); - - D = D + (alfa * residual) * u'; - - %uut=; - C = C - (alfa * u)* u'; - % lambda=(19*lambda+1)/20; - % DtD = DtD + alfa * ( v*u' + u*v') + alfa^2 * (residual'*residual) * uut; - -% if (mod(i,modi)==0) -% Ximd=zeros(size(X)); -% Ximd(:,p1((i-modi+1:i)))=data(:,i-modi+1:i); -% -% if (iternum==4) -% X_ima(:,:,1)=col2imstep(Ximd, [256 256], [8 8]); -% X_ima(:,:,2)=col2imstep(X, [256 256], [8 8]); -% X_ima(:,:,3)=zeros(256,256); -% else -% X_ima(:,:,1)=col2imstep(Ximd, [512 512], [8 8]); -% X_ima(:,:,2)=col2imstep(X, [512 512], [8 8]); -% X_ima(:,:,3)=zeros(512,512); -% end -% -% dictimg1=dictimg; -% dictimg = showdict(D,[8 8],... -% round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); -% dictimg1=(dictimg-dictimg1); -% -% figure(2); -% subplot(2,2,1); imshow(X_ima); title(sprintf('%d',i)); -% subplot(2,2,3); imshow(imresize(dictimg,2,'nearest')); -% subplot(2,2,4); imshow(imresize(dictimg1,2,'nearest')); -% subplot(2,2,2);imshow(C*(255/max(max(C)))); -% pause(0.02); -% if (i>=35000) -% modi=100; -% pause -% end; -% end -% end -end -%p1=p1(setxor(p2,1:end)); -%[D,cleared_atoms] = cleardict(D,X,muthresh,p1,replaced); -%replaced=zeros(dictsize,1); -% W=sparsecode(data, D, [], [], thresh); -% data=D*W; -%lambda=lambda+0.0002 -end -%Gamma=mexLasso(data, D, param); -%err=compute_err(D,Gamma, data); -%[y,i]=max(err); -%D(:,1)=data(:,i)/norm(data(:,i)); - -Dictionary = D;%D(:,p); -% figure(3); -% plot(lambda); -% mean(lambda); -% figure(4+j);plot(w_sp); -end - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% sparsecode % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -function Gamma = sparsecode(data,D,XtX,G,thresh) - -global CODE_SPARSITY codemode -global MEM_HIGH memusage -global ompfunc ompparams - -if (memusage < MEM_HIGH) - Gamma = ompfunc(D,data,G,thresh,ompparams{:}); - -else % memusage is high - - if (codemode == CODE_SPARSITY) - Gamma = ompfunc(D'*data,G,thresh,ompparams{:}); - - else - Gamma = ompfunc(D, data, G, thresh,ompparams{:}); - end - -end - -end - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% compute_err % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function err = compute_err(D,Gamma,data) - -global CODE_SPARSITY codemode - -if (codemode == CODE_SPARSITY) - err = sqrt(sum(reperror2(data,D,Gamma))/numel(data)); -else - err = nnz(Gamma)/size(data,2); -end - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% cleardict % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function [D,cleared_atoms] = cleardict(D,X,muthresh,unused_sigs,replaced_atoms) - -use_thresh = 4; % at least this number of samples must use the atom to be kept - -dictsize = size(D,2); - -% compute error in blocks to conserve memory -% err = zeros(1,size(X,2)); -% blocks = [1:3000:size(X,2) size(X,2)+1]; -% for i = 1:length(blocks)-1 -% err(blocks(i):blocks(i+1)-1) = sum((X(:,blocks(i):blocks(i+1)-1)-D*Gamma(:,blocks(i):blocks(i+1)-1)).^2); -% end - -cleared_atoms = 0; -usecount = replaced_atoms;%sum(abs(Gamma)>1e-7, 2); - -for j = 1:dictsize - - % compute G(:,j) - Gj = D'*D(:,j); - Gj(j) = 0; - - % replace atom - if ( (max(Gj.^2)>muthresh^2 || usecount(j)<use_thresh) && ~replaced_atoms(j) ) -% [y,i] = max(err(unused_sigs)); - D(:,j) = X(:,unused_sigs(end)) / norm(X(:,unused_sigs(end))); - unused_sigs = unused_sigs([1:end-1]); - cleared_atoms = cleared_atoms+1; - end -end - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% misc functions % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function err2 = reperror2(X,D,Gamma) - -% compute in blocks to conserve memory -err2 = zeros(1,size(X,2)); -blocksize = 2000; -for i = 1:blocksize:size(X,2) - blockids = i : min(i+blocksize-1,size(X,2)); - err2(blockids) = sum((X(:,blockids) - D*Gamma(:,blockids)).^2); -end - -end - - -function Y = colnorms_squared(X) - -% compute in blocks to conserve memory -Y = zeros(1,size(X,2)); -blocksize = 2000; -for i = 1:blocksize:size(X,2) - blockids = i : min(i+blocksize-1,size(X,2)); - Y(blockids) = sum(X(:,blockids).^2); -end - -end - -
--- a/DL/RLS-DLA/SMALL_rlsdla.m Tue Mar 15 15:30:54 2011 +0000 +++ b/DL/RLS-DLA/SMALL_rlsdla.m Wed Mar 16 13:41:02 2011 +0000 @@ -107,7 +107,7 @@ % Training data data=X; - +cnt=size(data,2); % C=(100000*thresh)*eye(dictsize); @@ -115,7 +115,7 @@ u=zeros(dictsize,1); -for i = 1:size(data,2) +for i = 1:cnt if (codemode == CODE_SPARSITY) w = ompmex(D,data(:,i),[],[],thresh,1,-1,0);
--- a/DL/RLS-DLA/SMALL_rlsdla1.m Tue Mar 15 15:30:54 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,479 +0,0 @@ -function Dictionary = SMALL_rlsdla1(X, params) - - - - - -global CODE_SPARSITY CODE_ERROR codemode -global MEM_LOW MEM_NORMAL MEM_HIGH memusage -global ompfunc ompparams exactsvd - -CODE_SPARSITY = 1; -CODE_ERROR = 2; - -MEM_LOW = 1; -MEM_NORMAL = 2; -MEM_HIGH = 3; - - -% p = randperm(size(X,2)); - - % coding mode % -X_norm=sqrt(sum(X.^2, 1)); -% X_norm_1=sum(abs(X)); -% X_norm_inf=max(abs(X)); -[X_norm_sort, p]=sort(X_norm);%, 'descend'); -% [X_norm_sort1, p5]=sort(X_norm_1);%, 'descend'); - -% if (isfield(params,'codemode')) -% switch lower(params.codemode) -% case 'sparsity' -% codemode = CODE_SPARSITY; -% thresh = params.Tdata; -% case 'error' -% codemode = CODE_ERROR; -% thresh = params.Edata; -% otherwise -% error('Invalid coding mode specified'); -% end -% elseif (isfield(params,'Tdata')) -% codemode = CODE_SPARSITY; -% thresh = params.Tdata; -% elseif (isfield(params,'Edata')) -% codemode = CODE_ERROR; -% thresh = params.Edata; -% -% else -% error('Data sparse-coding target not specified'); -% end - -thresh = params.Edata; -% max number of atoms % - -% if (codemode==CODE_ERROR && isfield(params,'maxatoms')) -% ompparams{end+1} = 'maxatoms'; -% ompparams{end+1} = params.maxatoms; -% end - - -% memory usage % - -if (isfield(params,'memusage')) - switch lower(params.memusage) - case 'low' - memusage = MEM_LOW; - case 'normal' - memusage = MEM_NORMAL; - case 'high' - memusage = MEM_HIGH; - otherwise - error('Invalid memory usage mode'); - end -else - memusage = MEM_NORMAL; -end - - -% iteration count % - -if (isfield(params,'iternum')) - iternum = params.iternum; -else - iternum = 10; -end - - -% omp function % - -if (codemode == CODE_SPARSITY) - ompfunc = @omp; -else - ompfunc = @omp2; -end - - -% % status messages % -% -% printiter = 0; -% printreplaced = 0; -% printerr = 0; -% printgerr = 0; -% -% verbose = 't'; -% msgdelta = -1; -% - -% -% for i = 1:length(verbose) -% switch lower(verbose(i)) -% case 'i' -% printiter = 1; -% case 'r' -% printiter = 1; -% printreplaced = 1; -% case 't' -% printiter = 1; -% printerr = 1; -% if (isfield(params,'testdata')) -% printgerr = 1; -% end -% end -% end -% -% if (msgdelta<=0 || isempty(verbose)) -% msgdelta = -1; -% end -% -% ompparams{end+1} = 'messages'; -% ompparams{end+1} = msgdelta; -% -% -% -% % compute error flag % -% -% comperr = (nargout>=3 || printerr); -% -% -% % validation flag % -% -% testgen = 0; -% if (isfield(params,'testdata')) -% testdata = params.testdata; -% if (nargout>=4 || printgerr) -% testgen = 1; -% end -% end - -% -% % data norms % -% -% XtX = []; XtXg = []; -% if (codemode==CODE_ERROR && memusage==MEM_HIGH) -% XtX = colnorms_squared(data); -% if (testgen) -% XtXg = colnorms_squared(testdata); -% end -% end - - -% mutual incoherence limit % - -if (isfield(params,'muthresh')) - muthresh = params.muthresh; -else - muthresh = 0.99; -end -if (muthresh < 0) - error('invalid muthresh value, must be non-negative'); -end - - - - - -% determine dictionary size % - -if (isfield(params,'initdict')) - if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:)))) - dictsize = length(params.initdict); - else - dictsize = size(params.initdict,2); - end -end -if (isfield(params,'dictsize')) % this superceedes the size determined by initdict - dictsize = params.dictsize; -end - -if (size(X,2) < dictsize) - error('Number of training signals is smaller than number of atoms to train'); -end - - -% initialize the dictionary % - -if (isfield(params,'initdict')) - if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:)))) - D1 = X(:,params.initdict(1:dictsize)); - D2 = X(:,params.initdict(1:dictsize)); - else - if (size(params.initdict,1)~=size(X,1) || size(params.initdict,2)<dictsize) - error('Invalid initial dictionary'); - end - D1 = params.initdict(:,1:dictsize); - D2 = params.initdict(:,1:dictsize); - end -else - data_ids = find(colnorms_squared(X) > 1e-6); % ensure no zero data elements are chosen - perm = randperm(length(data_ids)); - D = X(:,data_ids(perm(1:dictsize))); -end - - -% normalize the dictionary % - -% D = normcols(D); -% DtD=D'*D; - -err = zeros(1,iternum); -gerr = zeros(1,iternum); - -if (codemode == CODE_SPARSITY) - errstr = 'RMSE'; -else - errstr = 'mean atomnum'; -end -X(:,p(X_norm_sort<thresh))=0; -% if (iternum==4) -% X_im=col2imstep(X, [256 256], [8 8]); -% else -% X_im=col2imstep(X, [512 512], [8 8]); -% end -% figure(10); imshow(X_im); -p1=p(X_norm_sort>thresh); - -%p1=p1(p2(1:40000)); -%end-min(40000, end)+1:end));%1:min(40000, end))); -%p1 = randperm(size(data,2));%size(data,2) -%data=data(:,p1); - -C1=(100000*thresh)*eye(dictsize); -C2=(100000*thresh)*eye(dictsize); -% figure(11); -w=zeros(dictsize,1); -replaced=zeros(dictsize,1); -u=zeros(dictsize,1); - dictimg = showdict(D1,[8 8],round(sqrt(size(D1,2))),round(sqrt(size(D1,2))),'lines','highcontrast'); -% h=imshow(imresize(dictimg,2,'nearest')); -lambda=0.9998 -for j=1:3 -if size(p1,2)>20000 - p2 = randperm(floor(size(p1,2)/2)); - p2=sort(p2(1:20000)); - data1=X(:,p1(p2)); - data2=X(:,p1(floor(size(p1,2)/2)+p2)); -elseif size(p1,2)>0 - data=X(:,p1); -else - break; -end -% figure(1); -% plot(sqrt(sum(data.^2, 1))); -% a=size(data,2)/4; -% lambda0=0.99;%1-16/numS+iternum*0.0001-0.0002 -C1(1,1)=0; -C2(1,1)=0; -for i = 1:size(data1,2) -% if norm(data(:,i))>thresh - % par.multA= @(x,par) multMatr(D,x); % user function y=Ax - % par.multAt=@(x,par) multMatrAdj(D,x); % user function y=A'*x - % par.y=data(:,i); - % w=SolveFISTA(D,data(:,i),'lambda',0.5*thresh); - % w=sesoptn(zeros(dictsize,1),par.func_u, par.func_x, par.multA, par.multAt,options,par); - %w = SMALL_chol(D,data(:,i), 256,32, thresh);% - %w = sparsecode(data(:,i), D, [], [], thresh); - w1 = omp2mex(D1,data1(:,i),[],[],[],thresh,0,-1,-1,0); - w2 = omp2mex(D2,data2(:,i),[],[],[],thresh,0,-1,-1,0); - %w(find(w<1))=0; - %^2; -% lambda(i)=1-0.001/(1+i/a); -% if i<a -% lambda(i)=1-0.001*(1-(i/a)); -% else -% lambda(i)=1; -% end -% param.lambda=thresh; -% param.mode=2; -% param.L=32; -% w=mexLasso(data(:,i), D, param); - spind1=find(w1); - spind2=find(w2); - - %replaced(spind)=replaced(spind)+1; - %-0.001*(1/2)^(i/a); -% w_sp(i)=nnz(w); - residual1 = data1(:,i) - D1 * w1; - residual2 = data2(:,i) - D2 * w2; - %if ~isempty(spind) - %i - - C1 = C1 *(1/ lambda); - C2 = C2 *(1/ lambda); - u1 = C1(:,spind1) * w1(spind1); - u2 = C2(:,spind2) * w2(spind2); - %spindu=find(u); - % v = D' * residual; - - alfa1 = 1/(1 + w1' * u1); - alfa2 = 1/(1 + w2' * u2); - D1 = D1 + (alfa1 * residual1) * u1'; - D2 = D2 + (alfa2 * residual2) * u2'; - %uut=; - C1 = C1 - (alfa1 * u1)* u1'; - C2 = C2 - (alfa2 * u2)* u2'; - % lambda=(19*lambda+1)/20; - % DtD = DtD + alfa * ( v*u' + u*v') + alfa^2 * (residual'*residual) * uut; -% modi=5000; -% if (mod(i,modi)==0) -% Ximd=zeros(size(X)); -% Ximd(:,p((i-modi+1:i)))=data(:,i-modi+1:i); -% -% if (iternum==4) -% X_ima=col2imstep(Ximd, [256 256], [8 8]); -% else -% X_ima=col2imstep(Ximd, [512 512], [8 8]); -% end -% dictimg1=dictimg; -% dictimg = showdict(D,[8 8],... -% round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); -% dictimg1=(dictimg-dictimg1)*255; -% -% figure(2); -% subplot(2,2,1); imshow(X_ima); -% subplot(2,2,3); imshow(imresize(dictimg,2,'nearest')); -% subplot(2,2,4); imshow(imresize(dictimg1,2,'nearest')); -% subplot(2,2,2);imshow(C*(255/max(max(C)))); -% pause(0.02); -% end -% end -end -%p1=p1(setxor(p2,1:end)); -%[D,cleared_atoms] = cleardict(D,X,muthresh,p1,replaced); -%replaced=zeros(dictsize,1); -% W=sparsecode(data, D, [], [], thresh); -% data=D*W; -lambda=lambda+0.0001 -end -%Gamma=mexLasso(data, D, param); -%err=compute_err(D,Gamma, data); -%[y,i]=max(err); -%D(:,1)=data(:,i)/norm(data(:,i)); -% D=normcols(D); -% D_norm=sqrt(sum(D.^2, 1)); -% D_norm_1=sum(abs(D)); -% X_norm_1=sum(abs(X)); -% X_norm_inf=max(abs(X)); -% [D_norm_sort, p]=sort(D_norm_1, 'descend'); -Dictionary =[D1 D2]; -% figure(3); -% plot(lambda); -% mean(lambda); -% figure(4+j);plot(w_sp); -end - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% sparsecode % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -function Gamma = sparsecode(data,D,XtX,G,thresh) - -global CODE_SPARSITY codemode -global MEM_HIGH memusage -global ompfunc ompparams - -if (memusage < MEM_HIGH) - Gamma = ompfunc(D,data,G,thresh,ompparams{:}); - -else % memusage is high - - if (codemode == CODE_SPARSITY) - Gamma = ompfunc(D'*data,G,thresh,ompparams{:}); - - else - Gamma = ompfunc(D, data, G, thresh,ompparams{:}); - end - -end - -end - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% compute_err % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function err = compute_err(D,Gamma,data) - -global CODE_SPARSITY codemode - -if (codemode == CODE_SPARSITY) - err = sqrt(sum(reperror2(data,D,Gamma))/numel(data)); -else - err = nnz(Gamma)/size(data,2); -end - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% cleardict % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function [D,cleared_atoms] = cleardict(D,X,muthresh,unused_sigs,replaced_atoms) - -use_thresh = 4; % at least this number of samples must use the atom to be kept - -dictsize = size(D,2); - -% compute error in blocks to conserve memory -% err = zeros(1,size(X,2)); -% blocks = [1:3000:size(X,2) size(X,2)+1]; -% for i = 1:length(blocks)-1 -% err(blocks(i):blocks(i+1)-1) = sum((X(:,blocks(i):blocks(i+1)-1)-D*Gamma(:,blocks(i):blocks(i+1)-1)).^2); -% end - -cleared_atoms = 0; -usecount = replaced_atoms;%sum(abs(Gamma)>1e-7, 2); - -for j = 1:dictsize - - % compute G(:,j) - Gj = D'*D(:,j); - Gj(j) = 0; - - % replace atom - if ( (max(Gj.^2)>muthresh^2 || usecount(j)<use_thresh) && ~replaced_atoms(j) ) -% [y,i] = max(err(unused_sigs)); - D(:,j) = X(:,unused_sigs(end)) / norm(X(:,unused_sigs(end))); - unused_sigs = unused_sigs([1:end-1]); - cleared_atoms = cleared_atoms+1; - end -end - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% misc functions % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function err2 = reperror2(X,D,Gamma) - -% compute in blocks to conserve memory -err2 = zeros(1,size(X,2)); -blocksize = 2000; -for i = 1:blocksize:size(X,2) - blockids = i : min(i+blocksize-1,size(X,2)); - err2(blockids) = sum((X(:,blockids) - D*Gamma(:,blockids)).^2); -end - -end - - -function Y = colnorms_squared(X) - -% compute in blocks to conserve memory -Y = zeros(1,size(X,2)); -blocksize = 2000; -for i = 1:blocksize:size(X,2) - blockids = i : min(i+blocksize-1,size(X,2)); - Y(blockids) = sum(X(:,blockids).^2); -end - -end - -
--- a/DL/RLS-DLA/SMALL_rlsdlaFirstClustTry.m Tue Mar 15 15:30:54 2011 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,494 +0,0 @@ -function Dictionary = SMALL_rlsdla(X, params) - - - - - -global CODE_SPARSITY CODE_ERROR codemode -global MEM_LOW MEM_NORMAL MEM_HIGH memusage -global ompfunc ompparams exactsvd - -CODE_SPARSITY = 1; -CODE_ERROR = 2; - -MEM_LOW = 1; -MEM_NORMAL = 2; -MEM_HIGH = 3; - - -% p = randperm(size(X,2)); - - % coding mode % -X_norm=sqrt(sum(X.^2, 1)); -% X_norm_1=sum(abs(X)); -% X_norm_inf=max(abs(X)); -[X_norm_sort, p]=sort(X_norm);%, 'descend'); -% [X_norm_sort1, p5]=sort(X_norm_1);%, 'descend'); - -% if (isfield(params,'codemode')) -% switch lower(params.codemode) -% case 'sparsity' -% codemode = CODE_SPARSITY; -% thresh = params.Tdata; -% case 'error' -% codemode = CODE_ERROR; -% thresh = params.Edata; -% otherwise -% error('Invalid coding mode specified'); -% end -% elseif (isfield(params,'Tdata')) -% codemode = CODE_SPARSITY; -% thresh = params.Tdata; -% elseif (isfield(params,'Edata')) -% codemode = CODE_ERROR; -% thresh = params.Edata; -% -% else -% error('Data sparse-coding target not specified'); -% end - -thresh = params.Edata; -% max number of atoms % - -% if (codemode==CODE_ERROR && isfield(params,'maxatoms')) -% ompparams{end+1} = 'maxatoms'; -% ompparams{end+1} = params.maxatoms; -% end - - -% memory usage % - -if (isfield(params,'memusage')) - switch lower(params.memusage) - case 'low' - memusage = MEM_LOW; - case 'normal' - memusage = MEM_NORMAL; - case 'high' - memusage = MEM_HIGH; - otherwise - error('Invalid memory usage mode'); - end -else - memusage = MEM_NORMAL; -end - - -% iteration count % - -if (isfield(params,'iternum')) - iternum = params.iternum; -else - iternum = 10; -end - - -% omp function % - -if (codemode == CODE_SPARSITY) - ompfunc = @omp; -else - ompfunc = @omp2; -end - - -% % status messages % -% -% printiter = 0; -% printreplaced = 0; -% printerr = 0; -% printgerr = 0; -% -% verbose = 't'; -% msgdelta = -1; -% - -% -% for i = 1:length(verbose) -% switch lower(verbose(i)) -% case 'i' -% printiter = 1; -% case 'r' -% printiter = 1; -% printreplaced = 1; -% case 't' -% printiter = 1; -% printerr = 1; -% if (isfield(params,'testdata')) -% printgerr = 1; -% end -% end -% end -% -% if (msgdelta<=0 || isempty(verbose)) -% msgdelta = -1; -% end -% -% ompparams{end+1} = 'messages'; -% ompparams{end+1} = msgdelta; -% -% -% -% % compute error flag % -% -% comperr = (nargout>=3 || printerr); -% -% -% % validation flag % -% -% testgen = 0; -% if (isfield(params,'testdata')) -% testdata = params.testdata; -% if (nargout>=4 || printgerr) -% testgen = 1; -% end -% end - -% -% % data norms % -% -% XtX = []; XtXg = []; -% if (codemode==CODE_ERROR && memusage==MEM_HIGH) -% XtX = colnorms_squared(data); -% if (testgen) -% XtXg = colnorms_squared(testdata); -% end -% end - - -% mutual incoherence limit % - -if (isfield(params,'muthresh')) - muthresh = params.muthresh; -else - muthresh = 0.99; -end -if (muthresh < 0) - error('invalid muthresh value, must be non-negative'); -end - - - - - -% determine dictionary size % - -if (isfield(params,'initdict')) - if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:)))) - dictsize = length(params.initdict); - else - dictsize = size(params.initdict,2); - end -end -if (isfield(params,'dictsize')) % this superceedes the size determined by initdict - dictsize = params.dictsize; -end - -if (size(X,2) < dictsize) - error('Number of training signals is smaller than number of atoms to train'); -end - - -% initialize the dictionary % - -if (isfield(params,'initdict')) - if (any(size(params.initdict)==1) && all(iswhole(params.initdict(:)))) - D = X(:,params.initdict(1:dictsize)); - else - if (size(params.initdict,1)~=size(X,1) || size(params.initdict,2)<dictsize) - error('Invalid initial dictionary'); - end - D = params.initdict(:,1:dictsize); - end -else - data_ids = find(colnorms_squared(X) > 1e-6); % ensure no zero data elements are chosen - perm = randperm(length(data_ids)); - D = X(:,data_ids(perm(1:dictsize))); -end - -% normalize the dictionary % - -% D = normcols(D); -% DtD=D'*D; - -err = zeros(1,iternum); -gerr = zeros(1,iternum); - -if (codemode == CODE_SPARSITY) - errstr = 'RMSE'; -else - errstr = 'mean atomnum'; -end -%X(:,p(X_norm_sort<thresh))=0; -% if (iternum==4) -% X_im=col2imstep(X, [256 256], [8 8]); -% else -% X_im=col2imstep(X, [512 512], [8 8]); -% end -% figure(10); imshow(X_im); - -%p1=p(cumsum(X_norm_sort)./[1:size(X_norm_sort,2)]>thresh); -p1=p(X_norm_sort>thresh); -tic; idx=kmeans(X(:,p1)',4, 'Start', 'cluster','MaxIter',200); toc -D=[D D D D]; -dictsize1=4*dictsize; -% X(:,setxor(p1,1:end))=0; -% X_im=col2imstep(X, [256 256], [8 8]); -% figure(10); imshow(X_im); -% if iternum==2 -% D(:,1)=D(:,2); -% end -%p1=p1(p2(1:40000)); -%end-min(40000, end)+1:end));%1:min(40000, end))); -%p1 = randperm(size(data,2));%size(data,2) -%data=data(:,p1); - -C=(100000*thresh)*eye(dictsize1); -% figure(11); -w=zeros(dictsize,1); -replaced=zeros(dictsize,1); -u=zeros(dictsize,1); -% dictimg = showdict(D,[8 8],round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); -% figure(11);imshow(imresize(dictimg,2,'nearest')); -% pause(1); -lambda=0.99986;%3+0.0001*params.linc; -for j=1:1 -if size(p1,2)>60000 - p2 = randperm(size(p1,2)); - - p2=sort(p2(1:60000));%min(floor(size(p1,2)/2),40000))); - size(p2,2) - data=X(:,p1(p2)); -elseif size(p1,2)>0 - p2 = randperm(size(p1,2)); - size(p2,2) - data=X(:,p1); -else - break; -end -% figure(1); -% plot(sqrt(sum(data.^2, 1))); -% a=size(data,2)/4; -% lambda0=0.99;%1-16/numS+iternum*0.0001-0.0002 -%C(1,1)=0; -modi=1000; -for i = 1:size(data,2) -% if norm(data(:,i))>thresh - % par.multA= @(x,par) multMatr(D,x); % user function y=Ax - % par.multAt=@(x,par) multMatrAdj(D,x); % user function y=A'*x - % par.y=data(:,i); - % w=SolveFISTA(D,data(:,i),'lambda',0.5*thresh); - % w=sesoptn(zeros(dictsize,1),par.func_u, par.func_x, par.multA, par.multAt,options,par); - %w = SMALL_chol(D,data(:,i), 256,32, thresh);% - %w = sparsecode(data(:,i), D, [], [], thresh); - w = omp2mex(D(:,((idx(i)-1)*dictsize+1):idx(i)*dictsize),data(:,i),[],[],[],thresh,0,-1,-1,0); - - %w(find(w<1))=0; - %^2; -% lambda(i)=1-0.001/(1+i/a); -% if i<a -% lambda(i)=1-0.001*(1-(i/a)); -% else -% lambda(i)=1; -% end -% param.lambda=thresh; -% param.mode=2; -% param.L=32; -% w=mexLasso(data(:,i), D, param); - spind=find(w); - %replaced(spind)=replaced(spind)+1; - %-0.001*(1/2)^(i/a); -% w_sp(i)=nnz(w); - residual = data(:,i) - D (:,((idx(i)-1)*dictsize+1):idx(i)*dictsize)* w; - %if ~isempty(spind) - %i - if (j==1) - C = C *(1/ lambda); - end - u = C(((idx(i)-1)*dictsize+1):idx(i)*dictsize,((idx(i)-1)*dictsize)+spind) * w(spind); - - %spindu=find(u); - % v = D' * residual; - - alfa = 1/(1 + w' * u); - - D(:,((idx(i)-1)*dictsize+1):idx(i)*dictsize) = D (:,((idx(i)-1)*dictsize+1):idx(i)*dictsize)+ (alfa * residual) * u'; - - %uut=; - C (((idx(i)-1)*dictsize+1):idx(i)*dictsize,((idx(i)-1)*dictsize+1):idx(i)*dictsize)= C(((idx(i)-1)*dictsize+1):idx(i)*dictsize,((idx(i)-1)*dictsize+1):idx(i)*dictsize) - (alfa * u)* u'; - % lambda=(19*lambda+1)/20; - % DtD = DtD + alfa * ( v*u' + u*v') + alfa^2 * (residual'*residual) * uut; - -% if (mod(i,modi)==0) -% Ximd=zeros(size(X)); -% Ximd(:,p1((i-modi+1:i)))=data(:,i-modi+1:i); -% -% if (iternum==4) -% X_ima(:,:,1)=col2imstep(Ximd, [256 256], [8 8]); -% X_ima(:,:,2)=col2imstep(X, [256 256], [8 8]); -% X_ima(:,:,3)=zeros(256,256); -% else -% X_ima(:,:,1)=col2imstep(Ximd, [512 512], [8 8]); -% X_ima(:,:,2)=col2imstep(X, [512 512], [8 8]); -% X_ima(:,:,3)=zeros(512,512); -% end -% -% dictimg1=dictimg; -% dictimg = showdict(D,[8 8],... -% round(sqrt(size(D,2))),round(sqrt(size(D,2))),'lines','highcontrast'); -% dictimg1=(dictimg-dictimg1); -% -% figure(2); -% subplot(2,2,1); imshow(X_ima); title(sprintf('%d',i)); -% subplot(2,2,3); imshow(imresize(dictimg,2,'nearest')); -% subplot(2,2,4); imshow(imresize(dictimg1,2,'nearest')); -% subplot(2,2,2);imshow(C*(255/max(max(C)))); -% pause(0.02); -% if (i>=35000) -% modi=100; -% pause -% end; -% end -% end -end -%p1=p1(setxor(p2,1:end)); -%[D,cleared_atoms] = cleardict(D,X,muthresh,p1,replaced); -%replaced=zeros(dictsize,1); -% W=sparsecode(data, D, [], [], thresh); -% data=D*W; -lambda=lambda+0.0002 -end -%Gamma=mexLasso(data, D, param); -%err=compute_err(D,Gamma, data); -%[y,i]=max(err); -%D(:,1)=data(:,i)/norm(data(:,i)); -D=normcols(D); -D_norm=sqrt(sum(D.^2, 1)); -D_norm_1=sum(abs(D)); -% X_norm_1=sum(abs(X)); -% X_norm_inf=max(abs(X)); -[D_norm_sort, p]=sort(D_norm_1, 'descend'); -Dictionary = D;%D(:,p); -% figure(3); -% plot(lambda); -% mean(lambda); -% figure(4+j);plot(w_sp); -end - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% sparsecode % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -function Gamma = sparsecode(data,D,XtX,G,thresh) - -global CODE_SPARSITY codemode -global MEM_HIGH memusage -global ompfunc ompparams - -if (memusage < MEM_HIGH) - Gamma = ompfunc(D,data,G,thresh,ompparams{:}); - -else % memusage is high - - if (codemode == CODE_SPARSITY) - Gamma = ompfunc(D'*data,G,thresh,ompparams{:}); - - else - Gamma = ompfunc(D, data, G, thresh,ompparams{:}); - end - -end - -end - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% compute_err % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function err = compute_err(D,Gamma,data) - -global CODE_SPARSITY codemode - -if (codemode == CODE_SPARSITY) - err = sqrt(sum(reperror2(data,D,Gamma))/numel(data)); -else - err = nnz(Gamma)/size(data,2); -end - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% cleardict % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function [D,cleared_atoms] = cleardict(D,X,muthresh,unused_sigs,replaced_atoms) - -use_thresh = 4; % at least this number of samples must use the atom to be kept - -dictsize = size(D,2); - -% compute error in blocks to conserve memory -% err = zeros(1,size(X,2)); -% blocks = [1:3000:size(X,2) size(X,2)+1]; -% for i = 1:length(blocks)-1 -% err(blocks(i):blocks(i+1)-1) = sum((X(:,blocks(i):blocks(i+1)-1)-D*Gamma(:,blocks(i):blocks(i+1)-1)).^2); -% end - -cleared_atoms = 0; -usecount = replaced_atoms;%sum(abs(Gamma)>1e-7, 2); - -for j = 1:dictsize - - % compute G(:,j) - Gj = D'*D(:,j); - Gj(j) = 0; - - % replace atom - if ( (max(Gj.^2)>muthresh^2 || usecount(j)<use_thresh) && ~replaced_atoms(j) ) -% [y,i] = max(err(unused_sigs)); - D(:,j) = X(:,unused_sigs(end)) / norm(X(:,unused_sigs(end))); - unused_sigs = unused_sigs([1:end-1]); - cleared_atoms = cleared_atoms+1; - end -end - -end - - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% misc functions % -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - - -function err2 = reperror2(X,D,Gamma) - -% compute in blocks to conserve memory -err2 = zeros(1,size(X,2)); -blocksize = 2000; -for i = 1:blocksize:size(X,2) - blockids = i : min(i+blocksize-1,size(X,2)); - err2(blockids) = sum((X(:,blockids) - D*Gamma(:,blockids)).^2); -end - -end - - -function Y = colnorms_squared(X) - -% compute in blocks to conserve memory -Y = zeros(1,size(X,2)); -blocksize = 2000; -for i = 1:blocksize:size(X,2) - blockids = i : min(i+blocksize-1,size(X,2)); - Y(blockids) = sum(X(:,blockids).^2); -end - -end - -
--- a/Problems/generateImageDenoiseProblem.m Tue Mar 15 15:30:54 2011 +0000 +++ b/Problems/generateImageDenoiseProblem.m Wed Mar 16 13:41:02 2011 +0000 @@ -156,5 +156,5 @@ data.maxval = maxval; data.initdict= initdict; data.signalDim=2; - +data.sparse=1; end %% end of function \ No newline at end of file
--- a/To Do list.m Tue Mar 15 15:30:54 2011 +0000 +++ b/To Do list.m Wed Mar 16 13:41:02 2011 +0000 @@ -10,4 +10,12 @@ - adds noise sigma*var(signal) - do 2d circular fft per image and take random lines as specified by mask - mask is made to sample n*szt/fold lines in phase encode and time dimensions - - play movies of original, undersampled cardiac images and mask \ No newline at end of file + - play movies of original, undersampled cardiac images and mask + +SMALLboxSetup.m - Installation script: + +- Add make call for DL/RLS_DLA/private and Problems/private (still need to + shall these be private directories as in KSVD. Files inside are common + common routines that are used all around, so it makes sense to put them + in util directory) +- \ No newline at end of file
--- a/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m Tue Mar 15 15:30:54 2011 +0000 +++ b/examples/Image Denoising/SMALL_ImgDenoise_DL_test_KSVDvsRLSDLA.m Wed Mar 16 13:41:02 2011 +0000 @@ -42,7 +42,7 @@ % Problem. As an input we set the number of training patches. for noise_ind=1:1 for im_num=4:4 -SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',512, noise_level(noise_ind)); +SMALL.Problem = generateImageDenoiseProblem(test_image(im_num).i, 40000, '',256, noise_level(noise_ind)); SMALL.Problem.name=im_num; results(noise_ind,im_num).noisy_psnr=SMALL.Problem.noisy_psnr;