% alternates.m
%
% demonstrate some alternate matlab formulations of the regular and 
% irregular sampling and reconstruction matrices for the no aperture case
% to compare Dirichlet and FFT filtering-based schemes.  a number of
% different but equivalent ways to compute various matrices are included.
%
% written by D.G. Long at BYU 15 Nov 2015
%
% note: code is designed to illustrate the results, and not necessarily
% to be efficient. 
clear

% declare anonymous function
check_diff=@(S,A,B,s) fprintf('Check: %24s  [norm: %f]  max diff: %g\n',S,norm(A)*s,max(abs(A(:)-B(:))));

% (arbitrarily) set parameters for demo
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 samples

% compute d values
d1=N1/R1;  % must be an integer greater than one for these examples
d2=N2/R2;  % must be an integer greater than one for these examples

% convenient constants
N=N1*N2;
M=R1*R2;

% generate regular sampling locations using subroutine
slocu=randsample2d(Ru,N1,N2,M1,M2,1);  % array of indexes to regular samples

% directly compute regular sampling locations
uloc1=([1:R1]-1)*d1+1; % cols
uloc2=([1:R2]-1)*d2+1; % rows
[uind1,uind2]=meshgrid(uloc1,uloc2);
slocu1=sub2ind([N2 N1],uind2(:),uind1(:));

% alternate scheme for generating regular sample locations
slocu2=kron(ones(size(uloc1')),uloc2')+kron((uloc1'-1)*N2,ones(size(uloc2')));

% generate sample location matrix
Su=zeros([N,M]);
for k=1:M
  Su(slocu(k),k)=1;
end

% display locations of samples
myfigure(1)
imagesc(reshape(sum(Su,2),N2,N1))
title('Sample locations')
xlabel('N_1')
ylabel('N_2')

myfigure(2)
imagesc(Su)
title('Sample location matrix')
xlabel('sample number')
ylabel('pixel index')

myfigure(3); clf
for k=1:Ru
  ind=find(Su(:,k)==1);
  subplot(R2,R1,k)
  imagesc(reshape(Su(:,k),N2,N1))
  axis off
end

% create single 2d Dirichlet matrix using Dm2d subroutine
[x,y]=meshgrid((1:N1)-1,(1:N2)-1); % x=col, y=row
D2d=Dm2d(x,y,M1,M2,N1,N2);

% create single 2d Dirichlet matrix as Kronecker product using 
% ideal 1d lowpass filters
lpf1=zeros([N1,1]); lpf2=zeros([N2,1]);
lpf1(1:M1+1)=1; lpf1(N1-M1+1:N1)=1;
lpf2(1:M2+1)=1; lpf2(N2-M2+1:N2)=1;
% impulse at origin
D1=zeros(size(lpf1)); D1(1)=1;
D2=zeros(size(lpf2)); D2(1)=1;
% LPF'ed impulses (Dirichlet functions)
D1=ifft(lpf1.*fft(D1)); D2=ifft(lpf2.*fft(D2));
% generate matrix
D2d1=kron(D1',D2);
disp('note: max diff maybe non-zero due to numerical limitations')
check_diff('2d Dirichlet',D2d,D2d1,sqrt(d1*d2)) % check difference

% repeat using 1d Dirichlet function
D1=Dm(0:N1-1,M1,N1)';
D2=Dm(0:N2-1,M2,N2)';
D2d1=kron(D1',D2);
check_diff('2d Dirichlet',D2d,D2d1,sqrt(d1*d2)) % check difference

% create single 2d Dirichlet matrix as Kronecker product using 
% ideal 2d lowpass filtering
lpf=zeros([N2 N1]); 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;
% impulse at origin
tmp=zeros([N2 N1]); tmp(1,1)=1;
D2d2=ifft2(lpf.*fft2(tmp));
check_diff('2d Dirichlet again',D2d,D2d2,sqrt(d1*d2)) % check difference

% display 2d Dirichlet matrix
myfigure(4)
imagesc(D2d); colorbar
title('2D Dirichlet kernel at origin')
xlabel('N_1')
ylabel('N_2')

% display spectrum of 2d Dirichlet matrix
myfigure(5)
f1=0:N1-1; f2=0:N2-1;
F=fft2(D2d);
imagesc(f1,f2,abs(F)); colorbar
title(sprintf('Dirichlet kernel spectrum  %d(%d)x%d(%d)',N1,R1,N2,R2))
xlabel('Freq index (x 1/N_1)')
ylabel('Freq index (x 1/N_2)')

% generate convolution matrix using 1d shifts that converts samples 
% to fully sampled signal
C=zeros([N N]);
for n=1:N
  r=mod(n-1,N2)+1;
  c=floor((n-1)/N2)+1;
  D1a=shift(D1,c-1);  % 1d circulant shift of D1
  D2a=shift(D2,r-1);  % 1d circulant shift of D2
  tmp=kron(D1a',D2a);
  C(:,n)=tmp(:);
end

% show convolution matrix
myfigure(6)
imagesc(C);colorbar
title('Convolution matrix')

% generate convolution matrix using 2d shifts of 2d Dirichlet function
C1=zeros([N N]);
for n=1:N
  r=mod(n-1,N2)+1;
  c=floor((n-1)/N2)+1;
  tmp=shift2d(D2d,r-1,c-1);  % 2d circulant shift of D2d matrix
  C1(:,n)=tmp(:);
end
check_diff('Conv matrix check',C,C1,1) % check difference

% generate convolution matrix via kronecker product and 
% multiplication by matrices that circulantly shift 1d Dirichlet functions
D1=Dm(0:N1-1,M1,N1)';
D2=Dm(0:N2-1,M2,N2)';
C2=zeros([N N]);
for n=1:N
  c=floor((n-1)/N2)+1;
  Sh1=shift2d(eye([N1 N1]),c-1,0); % generate shift matrix
  D1a=Sh1*D1;  % shifted 1d Dirichlet function
  r=mod(n-1,N2)+1;
  Sh2=shift2d(eye([N2 N2]),r-1,0); % generate shift matrix
  D2a=Sh2*D2;  % shifted 1d Dirichlet function
  tmp=kron(D1a',D2a);
  C2(:,n)=tmp(:);
end
check_diff('Conv matrix check 2',C,C2,1) % check difference

% repeat using a 2d matrix multiplies
C2=zeros([N N]);
for n=1:N
  c=floor((n-1)/N2)+1;
  Sh1=shift2d(eye([N1 N1]),c-1,0); % generate shift matrix
  r=mod(n-1,N2)+1;
  Sh2=shift2d(eye([N2 N2]),r-1,0); % generate shift matrix
  ktmp=kron(Sh1',Sh2')';
  tmp=ktmp*D2d(:);
  C2(:,n)=tmp(:);
end
check_diff('Conv matrix check 3',C,C2,1) % check difference

% repeat using a single (huge) 2d matrix multiply
Q=zeros([N^2 N]);
for n=1:N
  c=floor((n-1)/N2)+1;
  Sh1=shift2d(eye([N1 N1]),c-1,0); % generate shift matrix
  r=mod(n-1,N2)+1;
  Sh2=shift2d(eye([N2 N2]),r-1,0); % generate shift matrix
  ktmp=kron(Sh1',Sh2')';
  Q((n-1)*N+1:n*N,:)=ktmp;
end
C2=reshape(Q*D2d(:),N,N);
check_diff('Conv matrix check 4',C,C2,1) % check difference

% directly generate convolution matrix using 2d Dirichlet kernal evaluations
C1=zeros([N N]);
for n=1:N
  r=mod(n-1,N2)+1;
  c=floor((n-1)/N2)+1;
  [x,y]=meshgrid(c-(1:N1),r-(1:N2));
  tmp=Dm2d(x,y,M1,M2,N1,N2);
  C1(:,n)=tmp(:);
end
check_diff('Conv matrix check again',C,C1,1) % check difference

% create Dirichlet-filtered sample locations via matrix multiplication
%
% note: since the sample locations are a regular cubic sampling, the
% result is D, the 2D regular reconstruction matrix
%
D=C*Su;

% create D using multiple 2d LPFing of sample locations
%  first, generate an ideal lowpass filter response
lpf=zeros([N2 N1]); 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;
%  then use FFTs to compute Dirichlet kernels
D1=zeros(size(Su));
for n=1:Ru
  tmp=reshape(Su(:,n),N2,N1);
  tmp2=ifft2(lpf.*fft2(tmp));
  D1(:,n)=tmp2(:);
end
check_diff('D matrix check',D,D1,sqrt(d1*d2)) % check difference

% directly create D using Dirichlet function evaluations
D2=zeros(size(Su));
for nn=1:Ru
  n=slocu(nn);
  r=mod(n-1,N2)+1;
  c=floor((n-1)/N2)+1;
  [x,y]=meshgrid((1:N1)-c,(1:N2)-r);
  tmp=Dm2d(x,y,M1,M2,N1,N2);
  D2(:,nn)=tmp(:);
end
check_diff('D matrix check again',D,D2,sqrt(d1*d2)) % check difference

% show D, the filtered sample locations
myfigure(7)
imagesc(D);colorbar
title('2d D matrix')
xlabel('sample index')
ylabel('pixel index')

% show rows of D
myfigure(8); clf
for k=1:Ru
  subplot(R2,R1,k)
  imagesc(reshape(D(:,k),N2,N1),[0.25/(d1*d2) 1/(d1*d2)])
  axis off
end

%norm(C-D*pinv(D)) % an interesting observation; pinv(D)*D=I
check_diff('C=D*pinv(D)',C,D*pinv(D),1) % check difference

%%%%%%%%%%%%%
%
% illustrate sampling and reconstruction using D
%

%
% first consider regular sampling and reconstruction
%
% create a random synthetic bandlimited image
True=ifft2(lpf.*fft2(rand([N2,N1])));

% (ideally) sample it
Usamples=True(slocu);

% reconstruct image from regular samples (inverse not needed since regular
% sampling is used)
URecon=reshape(D*Usamples,N2,N1)*d1*d2;

% compute worst-case reconsruction error
% (non-zero due to numerical limitations)
check_diff('Regular reconstruction',True,URecon,1/norm(True)) % check difference

%
% now illustrate irregular sampling and reconstruction
% first generate a full-rank random sampling
result=0;
while result==0  % make sure irregular sampling is full rank
  % generate random samples
  sloc=randsample2d(Ru,N1,N2,M1,M2,0);  % array of indexes to irrregular samples

  % check rank of reconstruction
  result=checkrankD(sloc,M1,M2,N1,N2);

  if result==0
    % disp(' (Random irregular sampling not full rank, trying again)')
  end
end

% create D_delta matrix from D matrix and locations of irregular samples
DD=D(sloc,:); % D_delta

% directly create D_delta using 1d Dirichlet function evaluations
[sloc2a,sloc1a]=ind2sub([N2 N1],sloc(:)); 
DD1=zeros([Ru,Ru]);
for nn=1:Ru
  DD1(:,nn)=Dm(uind1(nn)-sloc1a(:),M1,N1).*Dm(uind2(nn)-sloc2a(:),M2,N2);
end
check_diff('DD matrix check',DD,DD1,sqrt(d1*d2)) % check difference

% directly create D_delta using 2d Dirichlet function evaluations
[x1,x2]=meshgrid(sloc1a,uind1);
[y1,y2]=meshgrid(sloc2a,uind2);
DD2=Dm2d(x2-x1,y2-y1,M1,M2,N1,N2)';
check_diff('DD matrix check again',DD,DD2,sqrt(d1*d2)) % check difference

% create D_delta using convolution matrix C and irregular sample locations
S=zeros([N,M]); % irregular sample matrix
for k=1:M
  S(sloc(k),k)=1;
end
tmp=C*S;
DD3=tmp(slocu,:)'; % D_delta
check_diff('DD matrix check',DD,DD3,sqrt(d1*d2)) % check difference

% compute D_delta from convolution matrix C and sample matrices
DD4=S'*C*Su;
check_diff('DD matrix check',DD,DD4,sqrt(d1*d2)) % check difference
% this is equivalent to simply sampling C
tmp=C(slocu,:);
DD5=tmp(:,sloc)';
check_diff('DD matrix check',DD,DD5,sqrt(d1*d2)) % check difference

%
% irregularly sample the true image using ideal (delta function) window
samples=True(sloc);

% reconstruct image from samples
% this is done in two steps: first compute uniform samples using inv(DD)
% then multiply  by D to create fully-reconstructed function
tmp=inv(DD)*samples;
Recon=reshape(D*tmp,N2,N1);

% compute worst-case reconstruction error
% (non-zero due to numerical limitations)
check_diff('Irregular reconstruction',True,Recon,1/norm(True)) % check difference
check_diff('Reg. sample recon.',Usamples,tmp/(d1*d2),1/norm(Usamples)) % check difference
