% demo_irregrecon2dAp.m
%
% demonstrate discrete irregular reconstruction in two dimensions with
% variable apertures using irregrecon2dAp.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 30 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 bandlimit
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;

% 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 irregular sample locations
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


% select aperture cases for simulation
ApCase=7;         % aperture function test case (0..7, see below)
L1=2;             % aperture array support setting (size=2*L1+1)
L2=L1;            % aperture array support setting (size=2*L2+1)

% define apertures
Ap=zeros([2*L2+1,2*L1+1 Ri]);  % aperture matrix
Apshape=ones([1 Ri]); % keep track of which aperture used
Napshapes=1; % number of different apertures

switch ApCase
case 0  % fixed delta function
  for i=1:Ri
    Ap(L2+1,L1+1,i)=1;
  end
  ApCaseName='Fixed delta function';

case 1 % fixed rect window
  Ap(2:end-1,2:end-1,:)=1;
  ApCaseName='Fixed small rect window';

case 2 % fixed rect window with spectral null
       % note: this case has a null in the signal spectrum and so
       % cannot completely reconstruct the signal
  Ap(:,:,:)=1;
  ApCaseName='Fixed large rect window';

case 3  % fixed simple 5 term window
  for i=1:Ri
    Ap(L2+1,L1+1,i)=1;
    Ap(L2+1,L1,  i)=0.25;
    Ap(L2+1,L1+2,i)=0.25;
    Ap(L2,  L1+1,i)=0.25;
    Ap(L2+2,L1+1,i)=0.25;
  end
  ApCaseName='Fixed 2d 5 term window';

case 4 % fixed Fejer (triangle) window
  k1=1-abs(-L1:L1)/(L1+1);
  k2=1-abs(-L2:L2)/(L2+1);
  [kk1, kk2]=meshgrid(k1,k2);
  win=kk1.*kk2;
  for i=1:Ri
    Ap(:,:,i)=win+0.1; 
  end
  ApCaseName='Fixed 2d Fejer window';

case 5 % constant Hamming window case
  k1=pi*(-L1:L1)/L1;
  k2=pi*(-L2:L2)/L2;
  [kk1, kk2]=meshgrid(k1,k2);
  win=(0.56+0.48*cos(kk1)).*(0.56+0.48*cos(kk2));
  for i=1:Ri
    Ap(:,:,i)=win; 
  end
  Napshapes=1;
  ApCaseName='Constant 2d Hamming window';  

case 6 % variable aperture (randomly select two different Fejer windows)
  k1=1-abs(-L1:L1)/(L1+1);
  k2=1-abs(-L2:L2)/(L2+1);
  [kk1, kk2]=meshgrid(k1,k2);
  win1=kk1.*kk2;
  k1=[0 1-abs(1-L1:L1-1)/L1 0];
  k2=[0 1-abs(1-L2:L2-1)/L2 0];
  [kk1, kk2]=meshgrid(k1,k2);
  win2=kk1.*kk2;
  for i=1:Ri
    if rand>0.5
      Ap(:,:,i)=win1; 
      Apshape(i)=1;
    else
      Ap(:,:,i)=win2; 
      Apshape(i)=2;
    end
  end
  Napshapes=2;
  ApCaseName='Variable 2d Fejer windows';
  
case 7 % mix of Fejer and Hamming window cases
  k1=1-abs(-L1:L1)/(L1+1);
  k2=1-abs(-L2:L2)/(L2+1);
  [kk1, kk2]=meshgrid(k1,k2);
  win1=kk1.*kk2; 
  k1=[0 1-abs(1-L1:L1-1)/L1 0];
  k2=[0 1-abs(1-L2:L2-1)/L2 0];
  [kk1, kk2]=meshgrid(k1,k2);
  win2=kk1.*kk2;
  k1=pi*(-L1:L1)/L1;
  k2=pi*(-L2:L2)/L2;
  [kk1, kk2]=meshgrid(k1,k2);
  win3=(0.56+0.48*cos(kk1)).*(0.56+0.48*cos(kk2));
  for i=1:Ri
    rval=rand;
    if rval>0.35
      if rval>0.7
	Ap(:,:,i)=win1; 
	Apshape(i)=1;
      else
	Ap(:,:,i)=win2; 
	Apshape(i)=2;
      end
    else 
      Ap(:,:,i)=win3; 
      Apshape(i)=3;
    end
  end
  Napshapes=3;
  ApCaseName='Variable Fejer & Hamming windows';

otherwise % constant delta function
  for i=1:Ri
    Ap(L2+1,L1+1,i)=1;
  end
  ApCaseName='default constant delta function';
end
  
% normalize apertures
for i=1:Ri
  Ap(:,:,i)=Ap(:,:,i)/sum(sum(Ap(:,:,i)));
end

% display apertures
myfigure(7);clf
coltab='krbmcgb';
% sampling x,y coordinates
[sloc2a,sloc1a]=ind2sub([N2,N1],sloc); 
apmax=max(max(max(Ap)));
switch Napshapes
  case 1
    subplot(2,2,1)
    imagesc(squeeze(Ap(:,:,1)),[0 apmax]);colorbar
    axis image
    title(sprintf('Single aperture case: %s',ApCaseName));
    subplot(2,2,2)
    surf(squeeze(Ap(:,:,1)))
    caxis([0 apmax]); colorbar;
    subplot(2,2,3)
    ind=find(Apshape==1);
    if length(ind)>0
      Ap1=squeeze(Ap(:,:,ind(1)));
      tmp=zeros([N2,N1]);
      tmp(1:size(Ap1,1),1:1:size(Ap1,2))=Ap1;
      imagesc(abs(fft2(tmp)).^2)
      title('aperture spectrum');
      axis image
      colorbar;
    end
    subplot(2,2,4)
    imagesc(abs(fft2(Fy)).^2)
    title('signal spectrum');
    axis image
  case 2
    subplot(2,3,1)
    ind=find(Apshape==1);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]);colorbar
      axis image
      title(sprintf('Dual (1 %c) aperture case: %d',coltab(1),ApCase));
    end
    subplot(2,3,2)
    if length(ind)>0
      ind=find(Apshape==2);
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('Dual (2 %c) aperture case: %d',coltab(2),ApCase));
    end
    subplot(2,3,3)
    imagesc(zeros([N2,N1]),[-1 1])
    axis image
    hold on % show locations of irregular samples
    for k=1:length(sloc1a)
      plot(sloc1a(k),sloc2a(k),[coltab(Apshape(k)) '.']);
    end
    hold off
    %plot(Apshape,'*');title(ApCaseName);
    %xlabel('sample');ylabel('aperture number');
    subplot(2,3,4)
    ind=find(Apshape==1);
    if length(ind)>0
      Ap1=squeeze(Ap(:,:,ind(1)));
      tmp=zeros([N2,N1]);
      tmp(1:size(Ap1,1),1:1:size(Ap1,2))=Ap1;
      imagesc(abs(fft2(tmp)).^2)
      title('aperture spectrum');
      axis image
    end
    subplot(2,3,5)
    ind=find(Apshape==2);
    if length(ind)>0
      Ap1=squeeze(Ap(:,:,ind(1)));
      tmp=zeros([N2,N1]);
      tmp(1:size(Ap1,1),1:1:size(Ap1,2))=Ap1;
      imagesc(abs(fft2(tmp)).^2)
      title('aperture spectrum');
      axis image
    end
    subplot(2,3,6)
    imagesc(abs(fft2(Fy)).^2)
    title('signal spectrum');
    axis image
  case 3
    subplot(2,4,1)
    ind=find(Apshape==1);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]);colorbar
      axis image
      title(sprintf('Three (1 %c) aperture case: %d',coltab(1),ApCase));
    end
    subplot(2,4,2)
    ind=find(Apshape==2);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('Three (2 %c) aperture case: %d',coltab(2),ApCase));
    end
    subplot(2,4,3)
    ind=find(Apshape==3);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('Three (3 %c) aperture case: %d',coltab(3),ApCase));
    end
    subplot(2,4,4)
    imagesc(zeros([N2,N1]),[-1 1])
    axis image
    hold on % show locations of irregular samples
    for k=1:length(sloc1a)
      plot(sloc1a(k),sloc2a(k),[coltab(Apshape(k)) '.']);
    end
    hold off
    %plot(Apshape,'*');title(ApCaseName);
    %xlabel('sample');ylabel('aperture number');
    subplot(2,4,5)
    ind=find(Apshape==1);
    if length(ind)>0
      Ap1=squeeze(Ap(:,:,ind(1)));
      tmp=zeros([N2,N1]);
      tmp(1:size(Ap1,1),1:1:size(Ap1,2))=Ap1;
      imagesc(abs(fft2(tmp)).^2)
      title('aperture spectrum');
      axis image
    end
    subplot(2,4,6)
    ind=find(Apshape==2);
    if length(ind)>0
      Ap1=squeeze(Ap(:,:,ind(1)));
      tmp=zeros([N2,N1]);
      tmp(1:size(Ap1,1),1:1:size(Ap1,2))=Ap1;
      imagesc(abs(fft2(tmp)).^2)
      title('aperture spectrum');
      axis image
    end
    subplot(2,4,7)
    ind=find(Apshape==3);
    if length(ind)>0
      Ap1=squeeze(Ap(:,:,ind(1)));
      tmp=zeros([N2,N1]);
      tmp(1:size(Ap1,1),1:1:size(Ap1,2))=Ap1;
      imagesc(abs(fft2(tmp)).^2)
      title('aperture spectrum');
      axis image
    end
    subplot(2,4,8)
    imagesc(abs(fft2(Fy)).^2)
    title('signal spectrum');
    axis image
  
  
  case 4
    subplot(2,3,1)
    ind=find(Apshape==1);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]);colorbar
      axis image
      title(sprintf('Four (1 %c) aperture case: %d',coltab(1),ApCase));
    end
    subplot(2,3,2)
    ind=find(Apshape==2);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('Four (2 %c) aperture case: %d',coltab(2),ApCase));
    end
    subplot(2,3,3)
    ind=find(Apshape==3);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('Four (3 %c) aperture case: %d',coltab(3),ApCase));
    end
    subplot(2,3,4)
    ind=find(Apshape==4);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('Four (4 %c) aperture case: %d',coltab(4),ApCase));
    end
    subplot(2,3,5)
    imagesc(zeros([N2,N1]),[-1 1])
    axis image
    hold on % show locations of irregular samples
    for k=1:length(sloc1a)
      plot(sloc1a(k),sloc2a(k),[coltab(Apshape(k)) '.']);
    end
    hold off
    %plot(ApshapeY);title(ApCaseName);
    %xlabel('sample');ylabel('aperture number');
  otherwise
    subplot(2,3,1)
    ind=find(Apshape==1);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]);colorbar
      axis image
      title(sprintf('%d (1 %c) aperture case: %d',Napshapes,coltab(1),ApCase));
    end
    subplot(2,3,2)
    ind=find(Apshape==2);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('%d (2 %c) aperture case: %d',Napshapes,coltab(2),ApCase));
    end
    subplot(2,3,3)
    ind=find(Apshape==3);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;
      axis image
      title(sprintf('%d (3 %c) aperture case: %d',Napshapes,coltab(3),ApCase));
    end
    subplot(2,3,4)
    ind=find(Apshape==4);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;    
      axis image
      title(sprintf('%d (4 %c) aperture case: %d',Napshapes,coltab(4),ApCase));
    end
    subplot(2,3,5)
    ind=find(Apshape==5);
    if length(ind)>0
      imagesc(squeeze(Ap(:,:,ind(1))),[0 apmax]); colorbar;    
      axis image
      title(sprintf('%d (5 %c) aperture case: %d',Napshapes,coltab(5),ApCase));
    end
    subplot(2,3,6)
    imagesc(zeros([N2,N1]),[-1 1])
    axis image
    hold on % show locations of irregular samples
    for k=1:length(sloc1a)
      plot(sloc1a(k),sloc2a(k),[coltab(Apshape(k)) '.']);
    end
    hold off
    %plot(Apshape,'*');title(ApCaseName);
    %xlabel('sample');ylabel('aperture number');
end
%print -dpng fig7Ap.png

% test sampling w/apertures
[test2,Ka]=checkrankD_Ap(sloc,M1,M2,N1,N2,Ap);
if test2==1
  fprintf('Sampling realization with apertures is full-rank %d\n',Ka);
else
  fprintf('Sampling realization with apertures is NOT full-rank %d\n',Ka);
end

% sample original signal
fru=Fy(sloc);    % irregular samples w/o apertures
fs=Fy(slocu);    % regular samples (used for reference)

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

% generate aperture-filtered samples directly from original image
L=(size(Ap)-1)/2;
k1=-L(2):L(2); k2=-L(1):L(1);
[kk1, kk2]=meshgrid(k1,k2);
fr=zeros([Ri 1]);
for i=1:Ri
  X=mod(sloc1a(i)+kk1+N1-1,N1)+1;
  Y=mod(sloc2a(i)+kk2+N2-1,N2)+1;
  ind=sub2ind([N2 N1],Y,X);
  fr(i)=sum(sum(Fy(ind).*squeeze(Ap(:,:,i))));
  %squeeze(Ap(:,:,i))
  %[fru(i) fr(i) fru(i)-fr(i)]
end

% 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 irregular samples
plot(sloc1a,sloc2a,'k.');
hold off
axis image
colorbar;
title('Original signal with irregular 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('frequency index -N1/2:N1/2')
ylabel('frequency index -N2/2:N2/2')
title('magnitude squared Fourier Transform');
%print -dpng fig1Ap.png

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

%% call reconstruction routine w/o apertures for comparison
[frecon1 fu1 slocU1 D1u Du Ra K1]=irregrecon2d(fru,sloc,N1,N2,M1,M2);

% set reconstruction rank flag
if rank(D1)==Ru
  test0=1;
else
  test0=0;
end

myfigure(2);clf
% show original signal 
subplot(2,2,1)
Fymax=max(max(Fy)); Fymin=min(min(Fy));
imagesc(Fy,[Fymin Fymax])
%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');
axis image
colorbar;
subplot(2,2,2)
% show 2d surface
surf(Fy);caxis([Fymin Fymax]);
%imagesc(reshape(fru,R2,R1),[Fymin Fymax]);colorbar
title('Original signal');
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 % handle errors
  Frecon_max=Frecon_min+1;
end
imagesc(frecon,[Frecon_min Frecon_max])
title(sprintf('Irreg Recon w/aperture 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)
imagesc(frecon1,[Frecon_min Frecon_max])
title(sprintf('Irreg Recon w/no-aperture R=%d of %d',rank(D1u),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;
%print -dpng fig2Ap.png

myfigure(8); clf
subplot(2,2,1)
fs2=reshape(fs,R2,R1);
imagesc(fs2,[Frecon_min Frecon_max])
title('True regular sample');
axis image
colorbar;
subplot(2,2,2)
fu2=reshape(fu,R2,R1);
imagesc(fu2,[Frecon_min Frecon_max])
title('Recon w/aperture');
axis image
colorbar;
subplot(2,2,3)
fu2=reshape(fu1,R2,R1);
imagesc(fu2,[Frecon_min Frecon_max])
title('Recon w/no-aperture');
axis image
colorbar;
subplot(2,2,4)
imagesc(Fy,[Frecon_min Frecon_max])
title('True signal');
axis image
colorbar;
hold on % show locations of irregular samples
[sloc2U,sloc1U]=ind2sub([N2,N1],slocU); % lexicographic-ordered uniform sampling subs
plot(sloc1U,sloc2U,'k.');
hold off
%print -dpng fig8Ap.png

%
% compute and print summary errors
fprintf('Variable aperture irregular sampling demo\n');
fprintf(' Number of different apertures=%d  Case name: %s\n',Napshapes,ApCaseName);

fprintf('Maximum abs regular sample reconstruction error w/Ap:   %f\n',max(max(abs(fu-fs))));
fprintf('Maximum abs regular sample reconstruction error w/o Ap: %f\n',max(max(abs(fu1(:)-fs))));
fprintf('Maximum abs signal reconstruction error w/Ap:    %f\n',max(max(abs(frecon-Fy))));
fprintf('Maximum abs signal reconstruction error w/o Ap:  %f\n',max(max(abs(frecon1-Fy))));

RMS_recon=100*sqrt(sum((Fy-frecon).^2))/sqrt(sum(Fy.^2));
RMS_recon1=100*sqrt(sum((Fy-frecon1).^2))/sqrt(sum(Fy.^2));
fprintf('RMS signal reconst err in percent w/Ap:   %f  Rank %d of %d\n',RMS_recon,rank(D1),Ru);
fprintf('RMS signal reconst err in percent w/o Ap: %f  Rank %d of %d\n',RMS_recon1,rank(D1u),Ru);
fprintf('Ap Reconstruction matrix condition number:   D1=%f  Rank(D1)=%d of %d\n',cond(D1),rank(D1),Ru);
fprintf('noAp Reconstruction matrix condition number: D1=%f  Rank(D1)=%d of %d\n',cond(D1u),rank(D1u),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('No-aperture reconstruction has inadequate rank %d %d\n',Kr,Ru);

  % try recontruction again using regularization to improve estimates
  alph=1e-4; % arbitrary regularization parameter
  fprintf(' Trying Ap reconstruction with regularization parameter %f\n',alph);
  [frecon fu slocU D1 D Ra K]=irregrecon2dAp(fr,sloc,N1,N2,M1,M2,Ap,alph);

  myfigure(3)
  % 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)
  % show 2d surface
  surf(Fy);caxis([Fymin Fymax]);
  %colorbar
  title('Original signal');

  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   % handle errors
    Frecon_max=Frecon_min_1;
  end
  imagesc(frecon,[Frecon_min Frecon_max])
  title(sprintf('Regularized 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('regularized converted regular sample values');
  axis image
  colorbar;

  % compute and print summary errors
  fprintf(' Maximum abs regular sample reconstruction error w/Ap: %f\n',max(max(abs(fu-fs))));
  fprintf(' Maximum abs signal reconstruction error w/Ap:         %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\n',RMS_recon);
  
else
  myfigure(3); clf
  fprintf('Perfect reconstruction possible based on aperture reconstruction matrix rank %d %d\n',Kr,Ru);
end

% indicate if perfect reconstruction from sampling locations possible
% using the irregular-to-regular reconstruction matrix rank
Kr=round(rank(D1u));
if Kr<Ru
  fprintf('No-aperture reconstruction has inadequate rank %d %d\n',Kr,Ru);
  % try no-aperture recontruction again using regularization to improve estimates
  alph=1e-4; % arbitrary regularization parameter
  fprintf(' Trying no-aperture recontruction with regularization parameter %f\n',alph);
  [frecon1 fu1 slocU1 D1u Du Ra K1]=irregrecon2d(fru,sloc,N1,N2,M1,M2,alph);
  fprintf(' Maximum abs signal reconstruction error w/o Ap:  %f\n',max(max(abs(frecon1-Fy))));
  RMS_recon1=100*sqrt(sum((Fy-frecon1).^2))/sqrt(sum(Fy.^2));
  fprintf(' RMS signal reconst err in percent w/o Ap: %f  Rank %d of %d\n',RMS_recon1,rank(D1u),Ru);  
else
  fprintf('Perfect reconstruction possible based on no-aperture reconstruction matrix rank %d %d\n',Kr,Ru);
end

% optionally display reconstruction matrices
if 0
  myfigure(4);clf
  subplot(2,2,1)
  imagesc(D1);
  title(sprintf('D1 w/aperture Rank=%d',rank(D1)));
  axis image
  colorbar;
  
  subplot(2,2,2)
  imagesc(D1u);
  title(sprintf('D1 w/o aperture  Rank=%d',rank(D1u)));
  axis image
  colorbar;

  invD1=pinv(D1);
  subplot(2,2,3)
  imagesc(invD1);
  title(sprintf('inv D1 w/ap Cond=%d',round(cond(D1))));
  axis image
  colorbar;
  
  invD1u=pinv(D1u);
  subplot(2,2,4)
  imagesc(invD1u);
  title(sprintf('inv D1 w/o ap Cond=%d',round(cond(D1u))));
  axis image
  colorbar;
end