function sloc=randsample2d(R,N1,N2,M1,M2,rflag);
%
% function sloc=randsample2d(R,N1,N2,M1,M2,rflag)
%
% generate R randomly located samples in two dimension  R>=(2*M1+1)*(2*M2+1)
% ensuring that they are unique and within one period (1:N2,1:N1)
% sample locations are returned in lexicographic order.  To covvert
% to x,y location use [yloc,xloc]=ind2sub([N2 N1],sloc)
%
% if rflag==0, random locations with 2d random uniform distribution 
%              note: returns (2*M1+1)*(2*M2+1) sample locations, ignoring
%              the requested R number of sample locations
% if rflag==1, regularly spaced locations in 2d, ignores requested R
% if rflag==2, random locations with constrained 2d uniform distribution
% if rflag==3, locations gridded in 2d cubic lattice w/two 1d random
%              spacings. note: if length(R)==1,returns (2*M1+1)*(2*M2+1) sample 
%              locations, ignoring the requested R number of sample
%              locations; otherwise returns R(1)*R(2) locations
% if rflag==4, random locations in 2d curbit lattice w/two 1d constrained distri
%              note: if length(R)==1 returns (2*M1+1)*(2*M2+1) sample 
%              locations, ignoring the requested R number of sample
%              locations; otherwise returns R(1)*R(2) locations
%

% written by D.G. Long at BYU 22 Dec 2010
% revised by D.G. Long at BYU 01 Jan 2011 + added rflag=3,4 options

if (rflag==1 || rflag==3 || rflag==4) % cubic lattice gridded locations
  if (rflag==1) % uniformly spaced, gridded locations
                % overrides R, returns (2*M1+1)*(2*M2+1) locations
    loc1=([1:(2*M1+1)]-1)*(N1/(2*M1+1))+1;
    loc2=([1:(2*M2+1)]-1)*(N2/(2*M2+1))+1;
  else
    if rflag==3
      rflagt=0;
    else
      rflagt=2;
    end
    if length(R)==1  % overrides R, returns (2*M1+1)*2(M2+1) locations
      loc1=randsample1d(2*M1+1,N1,M1,rflagt);
      loc2=randsample1d(2*M2+1,N2,M2,rflagt);
    else             % uses specified 2d R
      loc1=randsample1d(R(1),N1,M1,rflagt);
      loc2=randsample1d(R(2),N2,M2,rflagt);
    end
  end
  [loc1,loc2]=meshgrid(loc1,loc2);
  sloc=sub2ind([N2 N1],loc2(:),loc1(:)); % results already sorted

else % random in 2d locations
    
  sloc=zeros([R 1]);
  it=floor([N2*rand N1*rand]);
  it(it<=0)=1;
  ir=1;
  sloc(ir)=sub2ind([N2 N1],it(1),it(2));
    
  if rflag==2  % constrained random case
    du=[N2/(2*M2+1) N1/(2*M1+1)]; % range of samples
    du(du<2)=2;
    it0=it;
    loc1=([1:(2*M1+1)]-1)*(N1/(2*M1+1))+1;
    loc2=([1:(2*M2+1)]-1)*(N2/(2*M2+1))+1;
    [loc1,loc2]=meshgrid(loc1,loc2);
    sloc2=sub2ind([N2 N1],loc2(:),loc1(:)); % uniform samples
    clear loc1 loc2
  end

  while ir<R
    iflag=1;
    icnt=0;
    lmul=1;
    while iflag>0
      if rflag==2 % constrained locations (uniform random distribution
                  % around uniformly spaced locations
        [it0(2) it0(1)]=ind2sub([N2,N1],sloc2(ir));
	it=it0+round(lmul*2*du.*(rand(size(du))-0.5));
      else        % independent, uniformly distributed random locations
	it=floor([N2*rand(1) N1*rand(1)]);
      end
      it(it<=0)=1; % ensure within range
      if it(1)>N2
	it(1)=N2;
      end
      if it(2)>N1
	it(2)=N1;
      end
      sloc1=sub2ind([N2 N1],it(1),it(2));
      in=find(sloc==sloc1); % ensure unique locations
      if isempty(in) % unique (disjoint) location
	ir=ir+1;
	sloc(ir)=sloc1;
	iflag=0;
      else           % duplicate location
	icnt=icnt+1;
      end
      if icnt>100    % helps ensure that disjoint locations can be generated
	lmul=lmul*2;
	icnt=0;
      end
    end
  end
  sloc=sort(sloc);   % sort sample location order for convenience 
end
