function p_zxy = kspacePlaneRecon(p_txy, dx, dy, dt, c, varargin)
%KSPACEPLANERECON   3D planar FFT reconstruction.
%
% DESCRIPTION:
%       kspacePlaneRecon takes an acoustic pressure time-series p_txy
%       recorded over an uniform array of sensor points on a plane, and
%       constructs an estimate of the initial acoustic pressure
%       distribution that gave rise to those measurements using an
%       algorithm based on the FFT. The input p_txy must be indexed
%       p_txy(time step, sensor x position, sensor y position), where the
%       sensor spacing is given by dx and dy, the temporal spacing given by
%       dt, and the sound speed in the propagation medium (which is assumed
%       to be acoustically homogeneous) is given by c. The output p_zxy is
%       indexed as p_zxy(sensor z position, sensor x position, sensor y
%       position).
% 
%       The code uses a k-space algorithm which performs (1) a Fourier
%       transform on the data p_tx along both t and x dimensions (into
%       wavenumber-frequency space, where the wavenumber component is along
%       the detector line), (2) a mapping, based on the dispersion relation
%       for a plane wave in an acoustically homogeneous medium, from
%       wavenumber-frequency space to wavenumber-wavenumber space (where
%       the second component is orthogonal to the detector line), and
%       finally (3) an inverse Fourier transform back from the wavenumber
%       domain to the spatial domain. The result is an estimate of the
%       initial acoustic pressure distribution from which the acoustic
%       waves originated. 
%
%       Steps (1) and (3) can be performed efficiently using the fast
%       Fourier transform (FFT); they are therefore fastest when the number
%       of samples and number of detector points are both powers of 2. The
%       mapping in step (2) requires an interpolation of the data from an
%       evenly spaced grid of points in the wavenumber-frequency domain to
%       an evenly-spaced grid of points in the wavenumber-wavenumber
%       domain. The option 'Interp' may be used to choose the interpolation
%       method.
%
%       The physics of photoacoustics requires that the acoustic pressure
%       is initially non-negative everywhere. The estimate of the initial
%       pressure distribution generated by this code may have negative
%       regions due to artefacts arising from differences between the
%       assumed model and the real situation, e.g., homogeneous medium vs.
%       real, somewhat heterogeneous, medium; infinite measurement surface
%       vs. finite-sized region-of-detection, etc. A positivity (or
%       non-negativity) condition can be enforced by setting the optional
%       'PosCond' to true which simply sets any negative parts of the final
%       image to zero. 
%
% USAGE:
%       p_zxy = kspacePlaneRecon(p_txy, dx, dy, dt, c)
%       p_zxy = kspacePlaneRecon(p_txy, dx, dy, dt, c, ...)
%
% INPUTS:
%       p_txy       - pressure time-series recorded over an evenly spaced
%                     array of sensor points on a line (indexed as x, t)
%       dx, dy      - spatial step [m]
%       dt          - time step [s]
%       c           - acoustically-homogeneous sound speed [m/s]
%
% OPTIONAL INPUTS:
%       Optional 'string', value pairs that may be used to modify the
%       default computational settings.
%
%       'Interp'    - String input controlling the interpolation method
%                     used by interp3 in the reconstruction (default =
%                     *nearest).
%       'PlotRecon' - Boolean controlling whether a plot of the
%                     reconstructed estimate of the initial acoustic
%                     pressure distribution is produced (default = false).
%       'PosCond'   - Boolean controlling whether a positivity condition is
%                     enforced on the reconstructed estimate of the initial
%                     acoustic pressure distribution (default = false).
%
% OUTPUTS:
%       p_zxy       - estimate of the initial acoustic pressure
%                     distribution (indexed as z, x, y)
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 2nd July 2009
%       last update - 17th July 2009
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
%
% See also interp3, kspaceLineRecon, makeGrid

% start timer
tic;

% define defaults
num_req_inputs = 5;
interp_method = '*nearest';
plot_recon = false;
positivity_cond = false;

% replace with user defined values if provided
if nargin < num_req_inputs
    error('Incorrect number of inputs');
elseif ~isempty(varargin)
    for input_index = 1:2:length(varargin)
        switch varargin{input_index}
            case 'Interp'
                interp_method = varargin{input_index + 1};            
            case 'PlotRecon'
                plot_recon = varargin{input_index + 1};
            case 'PosCond'
                positivity_cond = varargin{input_index + 1};                
            otherwise
                error('Unknown optional input');
        end
    end
end

% mirror the time domain data about t = 0 to allow the cosine transform to
% be computed using an FFT
p_txy = [flipdim(p_txy(2:end, :, :), 1); p_txy];

% extract the size of mirrored input data
[Nt, Nx, Ny] = size(p_txy);

% update command line status
disp('Running k-space planar reconstruction...');
disp(['  grid size: ' num2str(Nx) ' by ' num2str(Ny) ' by ' num2str((Nt+1)/2) ' pixels']);
disp(['  interpolation mode: ' interp_method]);

% create a computational grid that is evenly spaced in kx, ky and w, where 
% Nz = Nt and dz = dt*c
kgrid = makeGrid(Nx, dx, Ny, dy, Nt, dt*c);

% from the grid for kz, create a computational grid for w using the
% relation dz = dt*c; this represents the initial sampling of w 
w = c*kgrid.kz;

% create a second computational grid for w that is based on the real part
% of the dispersion relation w/c = (kx^2 + ky^2 + kz^2)^1/2. This gives an
% w grid that is evenly spaced in kz. This is used for the interpolation
% from p(w, kx, ky) to p(kz, kx, ky). Only real w is taken to force kz (and
% thus z) to be symmetrical about 0 after the interpolation.
w_new = (c*kgrid.k);

% compute the value of kz
kz = sqrt( (w/c).^2 - kgrid.kx.^2 );

% calculate the scaling factor manually replacing the DC value with its
% limit (otherwise NaN results)
sf = c^2*kz./(2*w);
sf(w == 0 & kgrid.kx == 0 & kgrid.ky == 0) = c/2;

% compute the FFT of the input data and scale
p_wkxky = sf.*fftshift(fftn(fftshift(p_txy)));

% exclude the inhomogeneous part of the wave
p_wkxky(abs(w) < abs(c*kgrid.kx)) = 0;

% compute the interpolation from p(w, kx, ky) to p(kz, kx, ky)
p_kzkxky = interp3(kgrid.kx, w, kgrid.ky, p_wkxky, kgrid.kx, w_new, kgrid.ky, interp_method);

% set values outside the interpolation range to zero
p_kzkxky(isnan(p_kzkxky)) = 0;

% compute the inverse FFT
p_zxy = real(ifftshift(ifftn(ifftshift(p_kzkxky))));

% remove the left part of the mirrored data which corresponds to the
% negative part of the mirrored time data
p_zxy = p_zxy( (Nt + 1)/2:Nt, :, :);

% correct the scaling - the forward FFT is computed with a spacing of dt
% and the reverse requires a spacing of dz = dt*c, the reconstruction
% assumes that p0 is symmetrical about z, and only half the plane collects
% data (first approximation to correcting the limited view problem)
p_zxy = 2*2*p_zxy./c;

% enfore positivity condition
if positivity_cond
    disp('  applying positivity condition...');
    p_zxy(p_zxy < 0) = 0;
end

% update command line status
disp(['  computation completed in ' scaleTime(toc)]);

% plot the reconstruction
if plot_recon
    
    % allocate axis dimensions
    x_axis = [0 Nx*dx];
    y_axis = [0 Ny*dy];
    z_axis = [0 (Nt/2)*dt*c];
    
    % select suitable axis scaling factor
    [x_sc scale prefix] = scaleSI(max([Nx*dx, Ny*dy, (Nt/2)*dt*c])); 
    
    % select suitable plot scaling factor
    plot_scale = max(p_zxy(:));
    
    % create the figures
    figure;
    subplot(2, 2, 1), imagesc(x_axis*scale, z_axis*scale, squeeze(p_zxy(:, :, round(end/2))), [-plot_scale, plot_scale]);
    axis image;
    title('x-z plane');
    subplot(2, 2, 2), imagesc(y_axis*scale, z_axis*scale, squeeze(p_zxy(:, round(end/2), :)), [-plot_scale, plot_scale]);
    axis image;
    xlabel(['(All axes in ' prefix 'm)']);
    title('y-z plane');
    subplot(2, 2, 3), imagesc(x_axis*scale, y_axis*scale, squeeze(p_zxy(round(end/2), :, :)), [-plot_scale, plot_scale]);
    axis image;
    title('x-y plane');
    colormap(getColorMap);
end