idamnjanovic@70: function y = sampgrid(x,blocksize,varargin)
idamnjanovic@70: %SAMPGRID Sample a multi-dimensional matrix on a regular grid.
idamnjanovic@70: %  Y = SAMPGRID(X,BLOCKSIZE,I1,I2,...,Ip) extracts block samples of size
idamnjanovic@70: %  BLOCKSIZE from the p-dimensional matrix X, arranging the samples as the
idamnjanovic@70: %  column vectors of the matrix Y. The locations of the (1,1,..,1)-th
idamnjanovic@70: %  elements of each block are given in the index vectors I1,I2,..Ip. The
idamnjanovic@70: %  total number of samples taken is length(I1)xlength(I2)x...xlength(Ip).
idamnjanovic@70: %  BLOCKSIZE should either be a p-element vector of the form [N1,N2,...Np],
idamnjanovic@70: %  or a scalar N which is shorthand for the square block size [N N ... N].
idamnjanovic@70: %
idamnjanovic@70: %  Example: Sample a set of blocks uniformly from a 2D image.
idamnjanovic@70: %
idamnjanovic@70: %    n = 512; blocknum = 20000; blocksize = [8 8];
idamnjanovic@70: %    im = rand(n,n);
idamnjanovic@70: %    [i1,i2] = reggrid(size(im)-blocksize+1, blocknum);
idamnjanovic@70: %    blocks = sampgrid(im, blocksize, i1, i2);
idamnjanovic@70: %
idamnjanovic@70: %  See also REGGRID.
idamnjanovic@70: 
idamnjanovic@70: %  Ron Rubinstein
idamnjanovic@70: %  Computer Science Department
idamnjanovic@70: %  Technion, Haifa 32000 Israel
idamnjanovic@70: %  ronrubin@cs
idamnjanovic@70: %
idamnjanovic@70: %  November 2007
idamnjanovic@70: 
idamnjanovic@70: 
idamnjanovic@70: p = ndims(x);
idamnjanovic@70: if (p==2 && any(size(x)==1) && length(blocksize)==1)
idamnjanovic@70:   p = 1;
idamnjanovic@70: end
idamnjanovic@70: 
idamnjanovic@70: if (numel(blocksize)==1)
idamnjanovic@70:   blocksize = ones(1,p)*blocksize;
idamnjanovic@70: end
idamnjanovic@70: 
idamnjanovic@70: n = zeros(1,p);
idamnjanovic@70: for i = 1:p
idamnjanovic@70:   n(i) = length(varargin{i});
idamnjanovic@70: end
idamnjanovic@70: 
idamnjanovic@70: nsamps = prod(n);
idamnjanovic@70: 
idamnjanovic@70: % create y of the same class as x
idamnjanovic@70: y = zeros(prod(blocksize),nsamps,class(x));
idamnjanovic@70: 
idamnjanovic@70: % ids() contains the index of the current block in I1..Ip
idamnjanovic@70: ids = ones(p,1);
idamnjanovic@70: 
idamnjanovic@70: % block_ids contains the indices of the current block in X
idamnjanovic@70: block_ids = cell(p,1);
idamnjanovic@70: for j = 1:p
idamnjanovic@70:   block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1;
idamnjanovic@70: end
idamnjanovic@70: 
idamnjanovic@70: for k = 1:nsamps
idamnjanovic@70:   block = x(block_ids{:});
idamnjanovic@70:   y(:,k) = block(:);
idamnjanovic@70:   
idamnjanovic@70:   % increment ids() and block_ids{}
idamnjanovic@70:   if (k<nsamps)
idamnjanovic@70:     j = 1;
idamnjanovic@70:     while (ids(j) == n(j))
idamnjanovic@70:       ids(j) = 1;
idamnjanovic@70:       block_ids{j} = varargin{j}(1) : varargin{j}(1)+blocksize(j)-1;
idamnjanovic@70:       j = j+1;
idamnjanovic@70:     end
idamnjanovic@70:     ids(j) = ids(j)+1;
idamnjanovic@70:     block_ids{j} = varargin{j}(ids(j)) : varargin{j}(ids(j))+blocksize(j)-1;
idamnjanovic@70:   end
idamnjanovic@70: end
idamnjanovic@70: