function [f varargout]=irregrecon2dAp(fr,xr,N1,N2,M1,M2,U,varargin)
%
% Illustrate reconstruction of a 2d periodic signal from irregular 
% aperture-filtered samples
%
%  Fy=irregrecon2dAp(fr,xr,N1,N2,M1,M2,U<,alpha>)
%
%  [Fy fu]=irregrecon2dAp...
%  [Fy fu xu]=irregrecon2dAp...
%  [Fy fu xu D1]=irregrecon2dAp...
%  [Fy fu xu D1 D]=irregrecon2dAp...
%  [Fy fu xu D1 D Ra]=irregrecon2dAp...
%  [Fy fu xu D1 D Ra K]=irregrecon2dAp...
%
% Discrete irregular 2d reconstruction (with variable apertures)
%
% inputs:
%  fr = vector of irregular samples of length R >= (2*M1+1)*(2*M2+1)
%  Xr = vector of locations of irregular samples where locations
%       are specified in lexicographic order
%  N1,N2 = Period (N1>=2*M1+1, N2>=2*M2+1)
%  M1,M2 = 2D signal bandlimit
%  U = aperture functions (L2,L1,R) where (L2<=N2,L1<=N1 odd) where
%      the apertures are regularly sampled at the high rate 
%  alpha = (optional) regularization parameter (def 0=no regularization)
%
% outputs:
%  Fy = array (N2,N1) of values of one period of reconstructed signal
%  fu = (optional) (2*M1+1)*(2*M2+1) samples of equivalent uniform samples
%  xu = (optional) (2*M1+1)*(2*M2+1) locations of equivalent uniform samples
%  D1 = (optional) matrix relating irregular samples to equiv uniform
%        samples fr=D1*fu 
%  D  = (optional) matrix relating equivent uniform samples to output
%        signal Fy=D*fu
%  Ra = (optional) rank of D1
%  K  = (optional) condition number of D1
%
% Basic theory:  Given a 2d periodic discrete signal Fy(n1,n2) with 
% period N1,N2 which is bandlimited to FFT2(Fy)(k1,k2)=0 for 
% |k1|>M1 or |k2|>M2, Fy(n1,n2) can be reconstructed using the formula
%
% Fy(n1,n2)=sum_{j=-M1}^M1 sum_{j=-M1}^M1 fu(j,k) d1 d2 
%                D_M1(n1 - j d1,N1) D_M2(n2 - k d2,N2)
%
% where d1=N1/(2M1+1), d2=N2/(2M2+1), and
% D_M1(n1,N1) = sin(pi (2M1+1) x/N1) / sin(pi x/N1) (Dirichlet function)
% D_M2(n2,N2) = sin(pi (2M2+1) x/N2) / sin(pi x/N2) (Dirichlet function)
%
% Given Ri (Ri>(2*M1+1)*(2*M2+1)) irregular samples of Fy(n1,n2) at
% disjoint locations xr (xr in (0...N1),(0...N2)) not all in a 2d 
% circulant line, then fu=pinv(D1)*fr and Fy=D*fu where D1 and D are 
% matrices of appropriately sampled Dirichlet functions.  
%
% When the samples (fu) represent aperture-filtered samples of the 
% signal, the theory can be extended even when the aperture functions
% used for each sample are different.  For aperture filtered samples,
% the entries in the D1 matrix are generated by convolving the aperture
% functions and the Dirichlet functions.
%
% Note: Unlike in the one-dimensional case where any locations will work, 
% some 2d samping configurations are not invertable. Even when they are,
% the recontruction matrix can be poorly conditions, particlarly for 
% large N and M, resulting in reconstruction errors.  Aperture filtered
% samples can improve both the rank and conditioning or degrade it
% depending on the sampling location and the aperture functions.
%
% see Long and Franz, Band-limited signal reconstruction from irregular
% samples with variable apertures, IEEE Trans. Geosci. Remote Sens., 2015
%

% written by D.G. Long at BYU 30 Dec 2010
% revised by D.G. Long at BYU 10 Sep 2011 + add regularization
% revised by D.G. Long at BYU 31 Dec 2014

if nargin > 7
  alph=varargin{1};
else
  alph=0;
end

if length(fr) ~= length(xr)
  error('Number of samples and their location must be the same length in irregrecon');
end

Ri=length(fr);      % number of irregular samples
Ru=(2*M1+1)*(2*M2+1); % number of equivalent uniform samples
d1=N1/(2*M1+1);     % equivalent uniform sample spacing, must be >= 1
d2=N2/(2*M2+1);     % equivalent uniform sample spacing, must be >= 1
if d1<1 || d2<1
  error(sprintf('d=N/(2*M+1) must be greater than one in irregrecon2d but N1=%d M1=%d and N2=%d M2=%d so d1=%f d2=%f',N1,M1,N2,M2,d1,d2));
end
if Ri<Ru
  error(sprintf('R=%d is less than (2*M1+1)*(2*M2+1) where M1=%d M2=%d, reconstruction impossible',R,M1,M2));
end

if size(U,1)>N2 || size(U,2)>N1
  error(sprintf('aperture array must be less than image dimensions'));
end
if size(U,3)~=Ri
  error(sprintf('number of apertures must be the same as number of irregular samples'));
end

% generate ideally spaced regular sampling locations
% note: this computation is exact only if d is an integer
sl1=[1:floor(d1):floor(d1)*(2*M1+1)];
sl2=[1:floor(d2):floor(d2)*(2*M2+1)];
[sl1a,sl2a]=meshgrid(sl1,sl2);
% lexicographic-ordered uniform sampling index
xu=sub2ind([N2 N1],sl2a,sl1a);   

% create uniform sampling reconstruction matrix
D=zeros([N1*N2 Ru]);
for c1=1:N1*N2
  [x2 x1]=ind2sub([N2 N1],c1);
  % following line equivalent to commented out code
  %for c2=1:Ru
  %  D(c1,c2)=d1*d2*Dm(x1-sl1a(c2),M1,N1)*Dm(x2-sl2a(c2),M2,N2);
  %end
  D(c1,:)=d1*d2*Dm(x1-sl1a(:)',M1,N1).*Dm(x2-sl2a(:)',M2,N2);
end

% x,y sample locations of irregular samples
[sloc2a,sloc1a]=ind2sub([N2 N1],xr); 

% generate reconstruction matrix D1 (w/o aperture)
% (converts irregular samples to uniform samples)
D1=zeros([Ri Ru]);
%%D1a=zeros([Ri Ru]);
for c1=1:Ru
  for c2=1:Ri
    % aperture case
    D1(c2,c1)=d1*d2*ApDm2d(sl1a(c1)-sloc1a(c2),sl2a(c1)-sloc2a(c2),M1,M2,N1,N2,squeeze(U(:,:,c2)));
    % no-aperture case using 2d Dirchlet function
    % D1(c2,c1)=d1*d2*Dm2d(sl1a(c1)-sloc1a(c2),sl2a(c1)-sloc2a(c2),M1,M2,N1,N2);
    % equivalent no-aperture case with 1-d Dirchlet function
    % D1(c2,c1)=d1*d2*Dm(sl1a(c1)-sloc1a(c2),M1,N1)*Dm(sl2a(c1)-sloc2a(c2),M2,N2);
  end
end

% convert irregular samples to regular samples and then reconstruct
% the signal from the regular samples

if alph==0
  % could use
  % f_vec=D*pinv(D1)*fr;
  % but rather than explicitly compute inverse of D1 
  % use matlab's linear system solver routine, which is more accurate
  fu=D1\fr;     % solve irreg linear system for equiv uniform samples
  f_vec=D*fu;   % signal reconstructed from regular samples
else
  % SVD-based solution formulation that includes Tikhonov regularization (slower)
  [U,S,V]=svd(D1);
  Sinv=zeros(size(S'));
  for k=1:min(size(S))
    Sinv(k,k)=S(k,k)/(S(k,k).^2+alph);
  end
  fu=V*Sinv*U'*fr; % solution to irreg linear system
  f_vec=D*fu;   % signal reconstructed from regular samples
end

% reshape uniform samples into matrix
%fu_matrix=reshape(fu,(2*M2+1),(2*M1+1));

% reshape result into output reconstructed 2d signal
f=reshape(f_vec,N2,N1); 

% output optional arguments
Nout=nargout-1;
if Nout>0
  varargout(1)={fu};
  if Nout>1
    varargout(2)={xu};
    if Nout>2
      varargout(3)={D1};
      if Nout>3
	varargout(4)={D};
	if Nout>4
	  varargout(5)={rank(D1)};
	  if Nout>5
	    varargout(6)={cond(D1)};
	  end
	end
      end
    end
  end
end
