vsolap2.m

Sparse Signal Representation using Overlapping Frames (matlab toolbox)

function W=VSolap2(X,Fin,S,VSalg,B,el)
% VSolap2    Vector Selection for overlapping frame and 2D signal
%            a larger frame is build by repeating the input frame several
% times, i.e actual frame is input frame repeated B times. 
% Then, for each corresponding signal block
% a number of vectors are selected using a block-oriented vector selection
% algorithm, VSxxx, given by VSalg (char array).
% The program works for overlapping frames and 2D signal.
% Remarks on size of variables:
% In this version N,K,P and B are given (by the size of) as input arguments
% All these numbers must be square numbers since
% the program use N1=N2=sqrt(N) P1=P2=sqrt(P) and B1=B2=sqrt(B)
% In a later version we could have X 4D, Fin 5D and [B1, B2] instead of B and
% [N1,N2,L1,L2]=size(X), [N1,N2,K,P1,P2]=size(Fin) to give N=N1*N2, P=P1*P2, B=B1*B2
%
% W=VSolap2(X,F,S,VSalg,B);       
% W=VSolap2(X,F,W,VSalg,B,el);       
%--------------------------------------------------------------------------------
% arguments:
%  W     - The weight matrix, W is a sparse matrix of size KxL = Kx(L1*L2)
%  X     - The signal reshaped into blocks of length N, size NxL = (N1*N2)x(L1*L2)
%  F     - The frame of synthesis vectors (dictionary), 
%          size NxKxP = (N1*N2)xKx(P1*P2) (P>1, and P a square number)
%          the vectors of F must be normalized, [F,nf]=NormalizeF(F);
%  S     - the third input argument may have different meaning depending on its size.
%    1x1   S is scalar, then it is assumed to be sparsness factor, 0<S<1,
%          a total of S*N*L weights will be selected by a GMP method.
%    1xL   number of vectors to select for each block of length N 
%          S(l) is number of non-zero weights in W(:,l)
%          Note: only when B==1, if B>1 the distribution of weights may be changed
%    KxL   and the third argument should be the previous weights, W.
%          Now S=full(sum(W~=0)); and used as previous case (size 1xL).
%          If VSalg='VSab2' previous weights are used as input when VSalg is called
%  VSalg - Which vector selection algorithm to use, it must be a function
%          called like 'w=VSfomp2(x,S);' and have F and FF as as global variables. 
%          Recommended (=tested) are 'VSfomp2', 'VSmp2' or 'VSab2'.
%  B     - number of blocks to use at each call to VSalg, the actual frame will
%          be a block-diagonal matrix of size N(B+P-1)xKB, where the input frame, 
%          reshaped into size NPxK, is used in the blocks on the diagonal.
%          We should have B larger than or equal to (P-1) to keep overlapping 
%          within the adjacent (large/expanded/actual) frame.
%  el    - extra loops to do, vector selection will be done (1+el) times 
%          for each extended block,   default is el=0
%--------------------------------------------------------------------------------
% Note: Arguments for this function is almost like for VSblock, and VSolap1

%--------------------------------------------------------------------------------
% Copyright (c) 2001.  Karl Skretting.  All rights reserved.
% Hogskolen in Stavanger (Stavanger University), Signal Processing Group
% Mail:  [e-mail]   Homepage:  http://www.ux.his.no/~karlsk/
% 
% HISTORY:  dd.mm.yyyy
% Ver. 1.0  04.10.2001  KS: made function based on BlockVS (based on FindW)
% Ver. 1.1  26.10.2001  KS: some minor changes
% Ver. 1.2  26.10.2001  KS: VSolap2 made based on VSolap1
%           05.11.2001  KS: VSolap2 seems to work (including the GMP part)
% Ver. 1.3  05.12.2002  KS: moved from ..\Frames2D\ to ..\FrameTools
%--------------------------------------------------------------------------------

global F FF 
Mfile='VSolap2';
Display=1;        % decide if function display information, 2 also waitbar

if (nargin<5)
   error([Mfile,' should have at least 5 arguments, see help.']);
end
[k,Ls]=size(S);
[n,L]=size(X);
[N,K,P]=size(Fin);           % note: input frame is called Fin, not F
if (N==64) & (P==4) & ((B==4) | (B==9)) & (L==4096)
   % N=64;K=128;P=4;B=4;L=4096;
   N1=sqrt(N);N2=sqrt(N);
   L1=sqrt(L);L2=sqrt(L);
   P1=sqrt(P);P2=sqrt(P);
   B1=sqrt(B);B2=sqrt(B);
else
   error([Mfile,' only ready for the special case ',...
         '((N==64) & (P==4) & (B==4,9) & (L==4096)).']);
end   
if P<2
   error([Mfile,' should have P>1, for P==1 use VSblock.']);
end
if (k==K) & (Ls==L)
   W=S;             % the old (previous) weights 
   S=full(sum(W~=0));
   [k,Ls]=size(S);
   WeightsGiven=1;
else
   W=sparse(K,L);    % the weight matrix is sparse and initial zeros
   WeightsGiven=0;
end
if k~=1
   error([Mfile,': size of third argument (S or W) is wrong, see help.']);
end
if n~=N
   error([Mfile,': size of X and F do not correspond, see help.']);
end   
if Ls==1   % S is sparseness factor, 0<S<1
   if 1
      % find the S*N*L non-zero weights by a GMP method
      disp('Savg given, Global Matching Pursuit (GMP) to find first weights.');
      % first define the indexes for shift of X and W (Q)
      % L=25;L1=5;L2=5;
      noshift=reshape(1:L,L1,L2);
      up=[noshift(2:L1,:);noshift(1,:)];
      left=[noshift(:,2:L2),noshift(:,1)];
      down=[noshift(L1,:);noshift(1:(L1-1),:)];
      right=[noshift(:,L2),noshift(:,1:(L2-1))];
      noshift=noshift(:)';   % make the indexes rows
      up=up(:)';down=down(:)';left=left(:)';right=right(:)';
      upleft=up(left);   % we may now combine these translations
      downright=down(right);
      % t=reshape(upleft,L1,L2)
      % [temp,I]=sort(upleft);  % gives I==downright
      % 
      R=X;F=full(Fin);
      SNL=S(1)*N*L;      % number of non-zeros in W left to select
      count=0;
      while SNL>0
         count=count+1;
         disp(['GMP: count=',int2str(count),'   SNL=',int2str(SNL)]);
         Q=zeros(K,L);               % without overlap: Q=F'*X
         Q=Q+F(:,:,1)'*R(:,noshift);
         Q=Q+F(:,:,2)'*R(:,up);
         Q=Q+F(:,:,3)'*R(:,left);
         Q=Q+F(:,:,4)'*R(:,upleft);
         %
         [QL,I]=max(abs(Q));
         [temp,II]=sort(QL);      % l=II(L) is the index for the largest
         for m=L:(-1):floor(24*L/25)
            l=II(m);k=I(l);       % Q(l,k) is a large weight
            if QL(l)              % must check special
               if ~W(k,l); SNL=SNL-1; end;
               if (SNL<0); break; end;        % terminate for-loop
               W(k,l)=W(k,l)+Q(k,l);          % update the weight
               % now update R
               R(:,noshift(l))=R(:,noshift(l))-F(:,k,1)*Q(k,l);
               R(:,up(l))=R(:,up(l))-F(:,k,2)*Q(k,l);
               R(:,left(l))=R(:,left(l))-F(:,k,3)*Q(k,l);
               R(:,upleft(l))=R(:,upleft(l))-F(:,k,4)*Q(k,l);
               % make sure the neighbors are not used in this loop
               neighbors=[up(l),right(l),down(l),left(l),up(right(l)),...
                     up(left(l)),down(right(l)),down(left(l)),noshift(l)];
               QL(neighbors)=0;
            end
         end
      end
      clear Q R noshift left right up down upleft downright SNL count temp
      disp('GMP ok');disp(' ');
      S=full(sum(W~=0));
      [k,Ls]=size(S);
      return
   else
      disp('Savg given, distribute weights evenly to find first weights.');
      temp=S;
      S=zeros(1,L);
      t1=0;
      for l=1:L
         t1=t1+N*temp;
         S(l)=floor(t1);
         t1=t1-S(l);
      end
      [k,Ls]=size(S);
   end
end
if Ls~=L
   error([Mfile,': size of X and S do not correspond, see help.']);
end   
if (nargout < 1); 
   error([Mfile,': function must have one output argument, see help.']); 
end

if (nargin<6); el=0; end;
if length(el)==0; el=0; end;
if length(B)==0; B=0; end;
if (el>25); el=25; end;   % not too many extra loops
if (B==0); B1=P1;B2=P2;B=B1*B2; end;  % 'minimum' size of frame

if (B<(P-1))
   disp([Mfile,': (B<(P-1)) is not good.']);
end   

% now input and output arguments are checked well, it should now be ok

% first Fin should be reshaped into size NxKxP1xP2
Fin=reshape(Fin,N,K,P1,P2);

% then, build the frame to use here
F=zeros(N*(B1+P1-1)*(B2+P2-1),K*B1*B2);
for b1=1:B1
   for b2=1:B2
      b=b1+(b2-1)*B1;
      for p1=0:(P1-1)
         for p2=0:(P2-1)
            bp=(b1+p1)+(b2+p2-1)*(B1+P1-1);
            F((1:N)+(bp-1)*N,(1:K)+(b-1)*K)=Fin(:,:,p1+1,p2+1);
         end
      end
   end
end
FF=F'*F;
%
Stot=sum(S);  % S is now 1xL
if Display
   disp([Mfile,': size F is ',int2str(size(F,1)),'x',int2str(size(F,2)),...
         ',  Savg=',int2str(floor(Stot*B/L)),...
         ',  L=',int2str(L),',  B=',int2str(B),',  P=',int2str(P),...
         ',  N=',int2str(N),',  K=',int2str(K),',  Stot=',int2str(Stot),'.']);
   hwbL = ceil((el+1)*L/B);
   t1=['Please wait while ',int2str(hwbL),' calls to ',VSalg,' is done.'];
   if Display>1
      hwbi = 0;
      hwb = waitbar(0,t1);
   else
      disp([Mfile,': ',t1]);
   end
   disp(' ');
end

xl=zeros(1,(B1+P1-1)*(B2+P2-1));    % the l-indexes in X
wl=zeros(1,B1*B2);                  % the l-indexes in W
for i_el=0:el
   temp=ceil(rand(1,1)*B1);         % random offset each time
   for l1=temp:B1:L1
      temp=ceil(rand(1,1)*B2);      % random offset each time
      for l2=temp:B2:L2
         % find wl and xl
         % also find xba,  the reconstructed signal for weights before and after      
         xba=zeros(N*(B1+P1-1)*(B2+P2-1),1); 
         for b1=0:(B1+P1-2)          % 0:2
            for b2=0:(B2+P2-2)       % 0:2 
               i1=l1+b1; if i1>L1; i1=i1-L1; end;     % 1:L1
               i2=l2+b2; if i2>L2; i2=i2-L2; end;     % 1:L2
               i=i1+(i2-1)*L1;       % 1:L
               if (b1<B1) & (b2<B2)
                  iw=b1+b2*B1+1;     % 1:4
                  wl(iw)=i;
               end
               ix=b1+b2*(B1+P1-1)+1; % 1:9
               xl(ix)=i;
               % here for xba
               xbai=(1:N)+(ix-1)*N;
               for p1=0:(P1-1)
                  for p2=0:(P2-1)
                     i1=l1+b1-p1; i1=mod(i1-1,L1)+1;      % 1:L1
                     i2=l2+b2-p2; i2=mod(i2-1,L2)+1;      % 1:L2
                     i=i1+(i2-1)*L1;                      % 1:L
                     if ~ismember(i,wl) 
                        xba(xbai)=xba(xbai)+Fin(:,:,p1+1,p2+1)*W(:,i);
                     end
                  end
               end
               %
            end
         end
         % xl, wl, and xba should be ok here
         w=W(:,wl);w=w(:);     
         x=X(:,xl);x=x(:);     
         s=sum(S(wl));         % number of weights available
         r=x-xba;
         %w0=w; r0=r-F*w;       % test/debug
         %
         if s==1
            % should only find one weight, then find the best
            w=zeros(size(w));
            c=(r'*F);                % the inner products
            [temp,i]=max(abs(c));i=i(1);
            w(i)=c(i);
         elseif s>1
            if strcmp(VSalg,'VSab2')
               if (full(sum(w~=0))==s)
                  w=VSab2(r,w);    %  Vector Selection using previous w
               else
                  w=VSfomp2(r,s);  %  Vector Selection selecting s vectors
                  % w=VSab2(r,w);    %  Vector Selection using previous w
               end
            elseif strcmp(VSalg,'VSfomp2')
               w=VSfomp2(r,s);  %  Vector Selection selecting s vectors
            else
               w=feval(VSalg,r,s);       % ! Vector Selection
            end
         else
            w=zeros(size(w));
         end
         %r1=r-F*w;if ((r0'*r0)<(r1'*r1)); disp(wl);w=w0; end;   % test/debug
         W(:,wl)=reshape(w,K,B);
         S(wl)=full(sum(W(:,wl)~=0));
         if sum(S(wl))<s
            % if fewer than available was used, we may select more in next block
            i1=l1;i2=l2+B2; 
            if (i2>L2); 
               i2=i2-L2; 
               i1=l1+B1;
               if (i1>L1); i1=i1-L1; end;
            end
            i=i1+(i2-1)*L1;
            S(i)=S(i)+s-sum(S(wl));
         end
         if Display>1
            hwbi = hwbi+1;
            waitbar(hwbi/hwbL,hwb);
         end
      end
   end
end
%
if Display>1
   close(hwb);
end

% check that S is correct, is this test relevant??
% temp=full(sum(W~=0));
% if sum(temp==S)~=L
%    disp([Mfile,': Logical error?, program did not track changes in S.']);
%    S=temp;
% end

Ssum=sum(S);
if Display
   disp([Mfile,': Number of weights selected is ',int2str(Ssum),'.']);
else
   if Ssum~=Stot
      disp([Mfile,': Number of weights used has changed from ',int2str(Stot),...
            ' to ',int2str(Ssum),'.']);
   end
end
%

return