idamnjanovic@1: function [varargout] = reggrid(sz,num,mode)
idamnjanovic@1: %REGGRID Regular sampling grid.
idamnjanovic@1: %  [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM) returns the indices
idamnjanovic@1: %  of a regular uniform sampling grid over a p-dimensional matrix with
idamnjanovic@1: %  dimensions N1xN2x...xNp. NUM is the minimal number of required samples,
idamnjanovic@1: %  and it is ensured that the actual number of samples, given by
idamnjanovic@1: %  length(I1)xlength(I2)x...xlength(Ip), is at least as large as NUM.
idamnjanovic@1: %
idamnjanovic@1: %  [I1,I2,...,Ip] = REGGRID([N1 N2 ... Np], NUM,'MODE') specifies the
idamnjanovic@1: %  method for distributing the samples along each dimension. Valid modes
idamnjanovic@1: %  include 'eqdist' (the default mode) and 'eqnum'. 'eqdist' indicates an
idamnjanovic@1: %  equal distance between the samples in each dimension, while 'eqnum'
idamnjanovic@1: %  indicates an equal number of samples in each dimension.
idamnjanovic@1: %
idamnjanovic@1: %  Notes about MODE:
idamnjanovic@1: %
idamnjanovic@1: %    1. The 'eqnum' mode will generally fail when the p-th root of NUM
idamnjanovic@1: %    (i.e. NUM^(1/p)) is larger than min([N1 N2 ... Np]). Thus 'eqdist' is
idamnjanovic@1: %    the more useful choice for sampling an arbitrary number of samples
idamnjanovic@1: %    from the matrix (up to the total number of matrix entries).
idamnjanovic@1: %  
idamnjanovic@1: %    2. In both modes, the equality (of the distance between samples, or
idamnjanovic@1: %    the number of samples in each dimension) is only approximate. This is
idamnjanovic@1: %    because REGGRID attempts to maintain the appropriate equality while at
idamnjanovic@1: %    the same time find a sampling pattern where the total number of
idamnjanovic@1: %    samples is as close as possible to NUM. In general, the larger {Ni}
idamnjanovic@1: %    and NUM are, the tighter the equality.
idamnjanovic@1: %
idamnjanovic@1: %  Example: Sample a set of blocks uniformly from a 2D image.
idamnjanovic@1: %
idamnjanovic@1: %    n = 512; blocknum = 20000; blocksize = [8 8];
idamnjanovic@1: %    im = rand(n,n);
idamnjanovic@1: %    [i1,i2] = reggrid(size(im)-blocksize+1, blocknum);
idamnjanovic@1: %    blocks = sampgrid(im, blocksize, i1, i2);
idamnjanovic@1: %
idamnjanovic@1: %  See also SAMPGRID.
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: %  November 2007
idamnjanovic@1: 
idamnjanovic@1: dim = length(sz);
idamnjanovic@1: 
idamnjanovic@1: if (nargin<3)
idamnjanovic@1:   mode = 'eqdist';
idamnjanovic@1: end
idamnjanovic@1: 
idamnjanovic@1: if (any(sz<1))
idamnjanovic@1:   error(['Invalid matrix size : [' num2str(sz) ']']);
idamnjanovic@1: end
idamnjanovic@1: 
idamnjanovic@1: if (num > prod(sz))
idamnjanovic@1:   warning(['Invalid number of samples, returning maximum number of samples.']);
idamnjanovic@1: elseif (num <= 0)
idamnjanovic@1:   if (num < 0)
idamnjanovic@1:     warning('Invalid number of samples, assuming 0 samples.');
idamnjanovic@1:   end
idamnjanovic@1:   for i = 1:length(sz)
idamnjanovic@1:     varargout{i} = [];
idamnjanovic@1:   end
idamnjanovic@1:   return;
idamnjanovic@1: end
idamnjanovic@1: 
idamnjanovic@1: 
idamnjanovic@1: if (strcmp(mode,'eqdist'))
idamnjanovic@1:   
idamnjanovic@1:   % approximate distance between samples: total volume divided by number of
idamnjanovic@1:   % samples gives the average volume per sample. then, taking the p-th root
idamnjanovic@1:   % gives the average distance between samples
idamnjanovic@1:   d = (prod(sz)/num)^(1/dim);
idamnjanovic@1:   
idamnjanovic@1:   % compute the initial guess for number of samples in each dimension.
idamnjanovic@1:   % then, while total number of samples is too large, decrese the number of
idamnjanovic@1:   % samples by one in the dimension where the samples are the most crowded.
idamnjanovic@1:   % finally, do the opposite process until just passing num, so the final
idamnjanovic@1:   % number of samples is the closest to num from above.
idamnjanovic@1:   
idamnjanovic@1:   n = min(max(round(sz/d),1),sz);   % set n so that it saturates at 1 and sz
idamnjanovic@1:   
idamnjanovic@1:   active_dims = find(n>1);    % dimensions where the sample num can be reduced
idamnjanovic@1:   while(prod(n)>num && ~isempty(active_dims))
idamnjanovic@1:     [y,id] = min((sz(active_dims)-1)./n(active_dims));
idamnjanovic@1:     n(active_dims(id)) = n(active_dims(id))-1;
idamnjanovic@1:     if (n(active_dims(id)) < 2)
idamnjanovic@1:       active_dims = find(n>1);
idamnjanovic@1:     end
idamnjanovic@1:   end
idamnjanovic@1: 
idamnjanovic@1:   active_dims = find(n<sz);    % dimensions where the sample num can be increased
idamnjanovic@1:   while(prod(n)<num && ~isempty(active_dims))
idamnjanovic@1:     [y,id] = max((sz(active_dims)-1)./n(active_dims));
idamnjanovic@1:     n(active_dims(id)) = n(active_dims(id))+1;
idamnjanovic@1:     if (n(active_dims(id)) >= sz(active_dims(id)))
idamnjanovic@1:       active_dims = find(n<sz);
idamnjanovic@1:     end
idamnjanovic@1:   end
idamnjanovic@1: 
idamnjanovic@1:   for i = 1:dim
idamnjanovic@1:     varargout{i} = round((1:n(i))/n(i)*sz(i));
idamnjanovic@1:     varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2);
idamnjanovic@1:   end
idamnjanovic@1:   
idamnjanovic@1: elseif (strcmp(mode,'eqnum'))
idamnjanovic@1:   
idamnjanovic@1:   % same idea as above
idamnjanovic@1:   n = min(max( ones(size(sz)) * round(num^(1/dim)) ,1),sz);
idamnjanovic@1: 
idamnjanovic@1:   active_dims = find(n>1);
idamnjanovic@1:   while(prod(n)>num && ~isempty(active_dims))
idamnjanovic@1:     [y,id] = min((sz(active_dims)-1)./n(active_dims));
idamnjanovic@1:     n(active_dims(id)) = n(active_dims(id))-1;
idamnjanovic@1:     if (n(active_dims(id)) < 2)
idamnjanovic@1:       active_dims = find(n>1);
idamnjanovic@1:     end
idamnjanovic@1:   end
idamnjanovic@1:   
idamnjanovic@1:   active_dims = find(n<sz);
idamnjanovic@1:   while(prod(n)<num && ~isempty(active_dims))
idamnjanovic@1:     [y,id] = max((sz(active_dims)-1)./n(active_dims));
idamnjanovic@1:     n(active_dims(id)) = n(active_dims(id))+1;
idamnjanovic@1:     if (n(active_dims(id)) >= sz(active_dims(id)))
idamnjanovic@1:       active_dims = find(n<sz);
idamnjanovic@1:     end
idamnjanovic@1:   end
idamnjanovic@1:   
idamnjanovic@1:   for i = 1:dim
idamnjanovic@1:     varargout{i} = round((1:n(i))/n(i)*sz(i));
idamnjanovic@1:     varargout{i} = varargout{i} - floor((varargout{i}(1)-1)/2);
idamnjanovic@1:   end
idamnjanovic@1: else
idamnjanovic@1:   error('Invalid sampling mode');
idamnjanovic@1: end
idamnjanovic@1: