function [result varargout]=checkrankD(sloc,M1,M2,N1,N2)
%
% result=checkrankD(sloc,M1,M2,N1,N2)
% [result,R]=checkrankD(sloc,M1,M2,N1,N2)
% [result,R,D1]=checkrankD(sloc,M1,M2,N1,N2)
%
% check sampling matrix rank for length(sloc)>=R1*R2 sample locations 
% store in sloc (sloc in 1..N1*N2) where R1=2*M1+1; R2=2*M2+1;
%
% inputs:
%  sloc  = array of unique sample indexes
%  N1,N2 = Period (N1>=2*M1+1, N2>=2*M2+1)
%  M1,M2 = 2D signal bandlimit
%
% outputs:
%  result = 1 if full-rank, =0 if not-full rank
%  rank   = <optional> actual matrix rank
%  D1     = <optional> constructed irregular to regular reconstruction matrix
%

% written by D.G. Long at BYU 29 Dec 2014

R1=2*M1+1;
R2=2*M2+1;
d1=N1/R1; % equivalent uniform sample spacing, must be >= 1, 
          % ideally should be an integer (see below)
d2=N2/R2; % equivalent uniform sample spacing, must be >= 1, 
          % ideally should be an integer (see below)
if d1<1 || d2<1
  error(sprintf('d=N/(2*M+1) must be greater than one but N1=%d M1=%d and N2=%d M2=%d so d1=%f d2=%f',N1,M1,N2,M2,d1,d2));
end
Ri=length(sloc);  % number of samples
if Ri<R1*R2
  %error(sprintf('R=%d is less than (2*M1+1)*(2*M2+1) where M1=%d M2=%d, exact reconstruction impossible',Ri,M1,M2));
  result=0;
  return;
end

% x,y sample locations of samples
[sloc2,sloc1]=ind2sub([N2 N1],sloc(:)); % equiv to following line
% sloc2=mod(sloc(:)-1,N2)+1; sloc1=floor((sloc(:)-1)/N2)+1;

% generate ideally-spaced regular sampling locations with first sample at 1,1
% note: this particular computation is exact only if d is an integer
sl1=[1:floor(d1):floor(d1)*R1];
sl2=[1:floor(d2):floor(d2)*R2];
[sl1a,sl2a]=meshgrid(sl1,sl2); % equiv to following line
%sl1a=repmat(sl1,R2,1); sl2a=repmat(sl2',1,R1);

% generate irregular sampling to regular sampling conversion matrix
D1=zeros([Ri R2*R1]);
for c1=1:R2*R1
  %following line is faster but is equivalent to:
  % for c2=1:Ri
  %   D1(c2,c1)=d1*d2*Dm(sl1a(c1)-sloc1a(c2),M1,N1)*Dm(sl2a(c1)-sloc2a(c2),M2,N2);
  % end
  D1(:,c1)=d1*d2*Dm(sl1a(c1)-sloc1(:),M1,N1).*Dm(sl2a(c1)-sloc2(:),M2,N2);
end

% compute sampling matrix rank
R=rank(D1);

% check for full-rank and set output flag
if R==R1*R2
  result=1;
else
  result=0;
end

if nargout > 1 % optionally return matrix rank
  varargout(1)={R};
  if nargout > 2 % optionally return matrix
    varargout(2)={D1};
  end
end
return