% demo_irregrecon2d.m
%
% demonstrate discrete irregular reconstruction in two dimensions 
% with ideal apertures using irregrecon2d.m
%
% includes options for different combinations of sampling and apertures
%
% generates a new random sampling for each run
%
% written by D.G. Long at BYU 29 Dec 2010 
% modifications by D.G. Long at BYU 25 Jan 2014
%
clear

% set parameters
M1=1; M2=2;               % low pass coefficients
R1=(2*M1+1);              % dimension of andlimit
R2=(2*M2+1);              % dimension of bandlimit
d_ideal1=3; d_ideal2=2;   % set perfect values of d
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
Ri=R1*R2;        % set number of irregular samples (Ri>=Ru)
% experimentally increase the number of irregular samples
%Ri=Ri+1;

% compute exact d values (not needed)
%d1=N1/R1;  % must be greater than one
%d2=N2/R2;  % must be greater than one

% select one option for sampling
rflag=0; % pure random sampling with a 2d uniform distibution
%rflag=1; % regular sampling
%rflag=2; % constrained random sampling
%rflag=3; % gridded cubic lattice random sampling
%rflag=4; % gridded cubic lattice constrained random sampling

% generate the sample and frequency indices
x1=[0:N1-1];  xf1=2*pi*x1/N1-pi;  xfn=x1-floor(N1/2);
x2=[0:N2-1];  xf2=2*pi*x2/N2-pi;  xfn=x2-floor(N2/2);
[xx1, xx2]=meshgrid(x1,x2);  % signal dimension (N2 x N1)

% generate sampling locations
sloc=randsample2d(Ri,N1,N2,M1,M2,rflag); % index to random irregular samples
slocu=randsample2d(Ru,N1,N2,M1,M2,1);    % index to regular samples for reference

% test sampling
[test1,K]=checkrankD(sloc,M1,M2,N1,N2);
if test1==1
  fprintf('Sampling realization is full-rank %d\n',K);
else
  fprintf('Sampling realization is NOT full-rank %d\n',K);
end

% generate an ideal lowpass filter response
lpf=zeros(size(xx1)); lpf(1:M2+1,1:M1+1)=1; lpf(N2-M2+1:N2,N1-M1+1:N1)=1;
lpf(1:M2+1,N1-M1+1:N1)=1; lpf(N2-M2+1:N2,1:M1+1)=1;

% generate a random signal and band limit it to (M2,M1)
% Fy is the "true" signal we want to reconstruct 
y=randn(size(xx1));          % random signal, not bandlimited
y=y+0.5*sin(2*pi*2*xx1/N1).*sin(2*pi*3*xx2/N2); % add sinusoid for fun
Fy=real(ifft2(fft2(y).*lpf));% apply filter to generate bandlimited signal

% generate sampled signal
fr=Fy(sloc);    % irregular samples
fs=Fy(slocu);   % regular samples (used for reference)

% sampling x,y coordinates
[sloc2a,sloc1a]=ind2sub([N2,N1],sloc); 

% plot the true signal and its spectrum
myfigure(1); clf
subplot(2,2,1)
surf(Fy)
title('Original signal')
subplot(2,2,2)
imagesc(Fy)
hold on % show locations of uniform samples
[sloc2U,sloc1U]=ind2sub([N2,N1],slocu); % lexicographic-ordered uniform sampling subs
plot(sloc1U,sloc2U,'k.');
hold off
axis image
colorbar;
title('Original signal with regular samples');
%title(sprintf('Period N=%d,%d  M=%d,%d',N1,N2,M1,M2));

subplot(2,2,3)
surf(abs(fft2(Fy)).^2)
title('Signal spectrum');
subplot(2,2,4)
imagesc(abs(fft2(Fy)).^2)
axis image
colorbar
xlabel('freq index -N1/2:N1/2')
ylabel('freq index -N2/2:N2/2')
title('Signal spectrum');

%print -dpng fig1.png

% call reconstruction routine with optional outputs
[frecon fu slocU D1 D K]=irregrecon2d(fr,sloc,N1,N2,M1,M2);

myfigure(2)
% show original signal 
subplot(2,2,1)
Fymax=max(max(Fy)); Fymin=min(min(Fy));
imagesc(Fy,[Fymin Fymax])
axis image
colorbar;
%title(sprintf('Period N=%d,%d  M=%d,%d',N1,N2,M1,M2));
hold on % show locations of irregular samples
plot(sloc1a,sloc2a,'k.');
hold off
title('Original signal with irregular sampling');

subplot(2,2,2)
fs1=reshape(fs,R2,R1);
imagesc(fs1,[Fymin Fymax])
axis image
colorbar;
title('Orignal regular samples');

subplot(2,2,3)
Frecon_min=min(min(frecon));
Frecon_max=max(max(frecon));
if isnan(Frecon_min)  % handle errors
  Frecon_min=-1;
end
if isnan(Frecon_max)  % handle errors
  Frecon_max=-1;
end
if Frecon_min == Frecon_max
  Frecon_max=Frecon_max+1;
end
imagesc(frecon,[Frecon_min Frecon_max])
title(sprintf('Irregular reconstructed  R=%d of %d',rank(D1),Ru));
axis image
hold on % show locations of equivalent uniform samples
[sloc2U,sloc1U]=ind2sub([N2,N1],slocU); % lexicographic-ordered uniform sampling subs
plot(sloc1U,sloc2U,'kx');
hold off
colorbar;

subplot(2,2,4)
fu1=reshape(fu,R2,R1);
imagesc(fu1,[Frecon_min Frecon_max])
title('Reconstructed regular samples');
axis image
colorbar;

%print -dpng fig2.png

% compute and print summary errors
fprintf('Maximum abs regular sample reconstruction error: %f\n',max(max(abs(fu-fs))));
fprintf('Maximum abs signal reconstruction error:         %f\n',max(max(abs(frecon-Fy))));

RMS_recon=100*sqrt(sum((Fy-frecon).^2))/sqrt(sum(Fy.^2));
fprintf('RMS signal reconstruction error in percent: %f  Rank %d of %d\n',RMS_recon,rank(D1),Ru);
fprintf('Reconstruction matrix condition number: D1=%f  Rank(D1)=%d of %d\n',cond(D1),rank(D1),Ru);

% indicate if perfect reconstruction from sampling locations possible
% using the irregular-to-regular reconstruction matrix rank
Kr=round(rank(D1));
if Kr<Ru
  fprintf('Perfect reconstruction NOT possible due to inadequate reconstruction matrix rank %d %d\n',Kr,Ru);
  test0=0;
  
  % try recontruction again using regularization to improve estimates
  alph=1e-4; % arbitrary regularization parameter
  fprintf(' Regularization parameter value used: %f\n',alph);
  [frecon fu slocU D1 D K]=irregrecon2d(fr,sloc,N1,N2,M1,M2,alph);

  myfigure(3);clf
  % show original signal 
  subplot(2,2,1)
  Fymax=max(max(Fy)); Fymin=min(min(Fy));
  imagesc(Fy,[Fymin Fymax])
  axis image
  colorbar;
  %title(sprintf('Period N=%d,%d  M=%d,%d',N1,N2,M1,M2));
  hold on % show locations of irregular samples
  plot(sloc1a,sloc2a,'k.');
  hold off
  title('Original signal with irregular sampling');

  subplot(2,2,2)
  fs1=reshape(fs,R2,R1);
  imagesc(fs1,[Fymin Fymax])
  title('Orig signal regular samples');
  axis image
  colorbar;
  
  subplot(2,2,3)
  Frecon_min=min(min(frecon));
  Frecon_max=max(max(frecon));
  if isnan(Frecon_min)  % handle errors
    Frecon_min=-1;
  end
  if isnan(Frecon_max)  % handle errors
    Frecon_max=-1;
  end
  imagesc(frecon,[Frecon_min Frecon_max])
  title(sprintf('Regularized Irreg recon  R=%d of %d',rank(D1),Ru));
  axis image
  hold on % show locations of equivalent uniform samples
  [sloc2U,sloc1U]=ind2sub([N2,N1],slocU); % lexicographic-ordered uniform sampling subs
  plot(sloc1U,sloc2U,'kx');
  hold off
  colorbar;

  subplot(2,2,4)
  fu1=reshape(fu,R2,R1);
  imagesc(fu1,[Frecon_min Frecon_max])
  title('Regularized recon reg samples');
  axis image
  colorbar;

  % compute and print summary errors
  fprintf(' Maximum abs regular sample reconstruction error: %f\n',max(max(abs(fu-fs))));
  fprintf(' Maximum abs signal reconstruction error:         %f\n',max(max(abs(frecon-Fy))));

  RMS_recon=100*sqrt(sum((Fy-frecon).^2))/sqrt(sum(Fy.^2));
  fprintf(' RMS signal reconstruction error in percent: %f  Rank %d of %d\n',RMS_recon,rank(D1),Ru);

else
  myfigure(3); clf
  fprintf('Perfect reconstruction is possible based on reconstruction matrix rank %d of %d\n',Kr,Ru);
  test0=1;
end


% optionally display reconstruction matrices
if 0
  myfigure(4);clf
  subplot(1,2,1)
  surf(D1)
  subplot(1,2,2)
  imagesc(D1);
  title(sprintf('D1: irregular to regular conv R=%d of %d',rank(D1),Ru));
  axis image
  colorbar;
  
  myfigure(5);clf
  invD1=pinv(D1);
  subplot(1,2,1)
  surf(invD1)
  subplot(1,2,2)
  imagesc(invD1);
  title(sprintf('inv D1  Cond=%d',round(cond(D1))));
  axis image
  colorbar;

  myfigure(6);clf
  subplot(1,2,1)
  surf(D)
  subplot(1,2,2)
  imagesc(D);
  title(sprintf('D for uniform sampling %d by %d',R1,R2));
  colorbar;
end