idamnjanovic@1: function cnt = countcover(sz,blocksize,stepsize)
idamnjanovic@1: %COUNTCOVER Covering of signal samples by blocks
idamnjanovic@1: %  CNT = COUNTCOVER(SZ,BLOCKSIZE,STEPSIZE) assumes a p-dimensional signal
idamnjanovic@1: %  of size SZ=[N1 N2 ... Np] covered by (possibly overlapping) blocks of
idamnjanovic@1: %  size BLOCKSIZE=[M1 M2 ... Mp]. The blocks start at position (1,1,..,1)
idamnjanovic@1: %  and are shifted between them by steps of size STEPSIZE=[S1 S2 ... Sp].
idamnjanovic@1: %  COUNTCOVER returns a matrix the same size as the signal, containing in
idamnjanovic@1: %  each entry the number of blocks covering that sample.
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: %  Ron Rubinstein
idamnjanovic@1: %  Computer Science Department
idamnjanovic@1: %  Technion, Haifa 32000 Israel
idamnjanovic@1: %  ronrubin@cs
idamnjanovic@1: %
idamnjanovic@1: %  August 2008
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: cnt = ones(sz);
idamnjanovic@1: for k = 1:length(sz)
idamnjanovic@1:   
idamnjanovic@1:   % this code is modified from function NDGRID, so it computes one
idamnjanovic@1:   % output argument of NDGRID at a time (to conserve memory)
idamnjanovic@1:   ids = (1:sz(k))';
idamnjanovic@1:   s = sz; s(k) = [];
idamnjanovic@1:   ids = reshape(ids(:,ones(1,prod(s))),[length(ids) s]);
idamnjanovic@1:   ids = permute(ids,[2:k 1 k+1:length(sz)]);
idamnjanovic@1:   
idamnjanovic@1:   cnt = cnt .* max( min(floor((ids-1)/stepsize(k)),floor((sz(k)-blocksize(k))/stepsize(k))) - ...
idamnjanovic@1:                     max(ceil((ids-blocksize(k))/stepsize(k)),0) + 1 , 0 );
idamnjanovic@1: end