%
% test code for evaluating sampling invertability for 2d sampling w/no aperture
%
% given a set of parameters (set below), generates a random or regular 
% 2d sampling and then tests the sampling to see if it is fully
% reconstructable.  Optionally plots the sampling and errors.  
% rerun for a new realization of the random sampling
%


%
% written by D.G. Long at BYU 29 Dec 2010 
%
clear

PLOT=1;  % include plots if set to 1

% set parameters
M1=2; M2=3;               % low pass coefficients
R1=(2*M1+1);              % dimension of bandlimit
R2=(2*M2+1);              % dimension of bandlimit
d_ideal1=3; d_ideal2=2;   % set perfect values of d

Ri=R1*R2;  % set number of irregular samples (Ri>=Ru)
%Ri=Ri+1;  % optionally increase the number of irregular samples

N1=R1*d_ideal1;  % signal period (for exact uniform spacing)
N2=R2*d_ideal2;  % signal period (for exact uniform spacing)
Ru=R1*R2;        % ideal number of regular samples

% select one option for sampling
rflag=0; % pure random sampling
%rflag=1; % regular sampling
%rflag=2; % constrained random sampling
%rflag=3; % gridded random sampling
%rflag=4; % gridded constrained random sampling

%%%%%%%%%%%%%

% generate sampling locations
sloc=randsample2d(Ri,N1,N2,M1,M2,rflag); % index to random irregular samples

% convert list to x,y indices
[sloc1,sloc2]=ind2sub([N1 N2],sloc(:)); % equiv to following line
sloc1a=mod(sloc(:)-1,N1)+1; sloc2a=floor((sloc(:)-1)/N1)+1;

% check to see if sampling is full-rank
test1=checkrankD(sloc,M1,M2,N1,N2);

% generate sampled function values for a constant surface (all 1's)
CONST=1;                    % true signal
fr=zeros(size(sloc))+CONST; % sample values

% call 2d signal reconstruction routine with optional outputs
% to return reconstructed signal and matrices
[frecon fu slocU D1 D K]=irregrecon2d(fr,sloc,N1,N2,M1,M2,0);

if PLOT % optionally plot sample locations
  myfigure(1)
  plotlocs(sloc,M1,M2,N1,N2);
  if test1
    title(sprintf('Full rank sampling  Cond=%f',cond(D1)))
  else
    title('Sampling not full rank')
  end
end

if PLOT % optionally plot the reconstruction error
  myfigure(2)
  imagesc(frecon-CONST);colorbar
  [sloc2 sloc1]=ind2sub([N2 N1],sloc);
  hold on;
  plot(sloc1,sloc2,'k.');
  hold off;
  title(sprintf('Recon Error mean=%f std=%f',mean(frecon(:))-1,std(frecon(:)-1)));
end

% indicate if perfect reconstruction from sampling locations possible
% using the irregular-to-regular reconstruction matrix rank
Kr=round(rank(D1));
fprintf('Reconstruction matrix:\n  Condition number(D1)=%f\n  Rank(D1)=%d of %d\n',cond(D1),Kr,Ru);
fprintf('  Mean error: %f\n  STD error: %f\n',mean(frecon(:)-CONST),std(frecon(:)-CONST));
if Kr<Ru
  fprintf('  Perfect reconstruction NOT possible. Inadequate matrix rank %d %d\n',Kr,Ru);
  test0=0;
else
  fprintf('  Perfect reconstruction possible based on matrix rank %d %d\n',Kr,Ru);
  test0=1;
end

if ~test0 % if not reconstructable, call reconstruction routine with
          % regularization parameter
  alph=0.001; % arbitrary regularization parameter
  [frecon fu slocU D1]=irregrecon2d(fr,sloc,N1,N2,M1,M2,alph);
  fprintf('  Regularized reconstruction with alpha=%f\n',alph);
  fprintf('    Mean error: %f\n    STD error: %f\n',mean(frecon(:)-CONST),std(frecon(:)-CONST));

  if PLOT % optionally plot the reconstruction error
    myfigure(3)
    imagesc(frecon-1);colorbar
    [sloc2 sloc1]=ind2sub([N2 N1],sloc);
    hold on;
    plot(sloc1,sloc2,'k.');
    hold off;
    title(sprintf('Reg. Recon Error mean=%f std=%f',mean(frecon(:)-CONST),std(frecon(:)-CONST)));
  end
  
else
    if PLOT % optionally plot the reconstruction error
      myfigure(3);clf
    end
end

