function p = kspaceLineRecon(p, dy, dt, c, varargin)
%KSPACELINERECON    2D linear FFT reconstruction.
%
% DESCRIPTION:
%       kspaceLineRecon takes an acoustic pressure time-series p_ty
%       recorded over an evenly spaced array of sensor points on a line,
%       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_ty must be indexed
%       p_ty(time step, sensor y position), where the sensor spacing is given
%       by 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_xy is indexed as
%       p_xy(x position, y position).
%
%       The code uses a k-space algorithm which performs (1) a Fourier
%       transform on the data p_ty along both t and y 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_xy = kspaceLineRecon(p_ty, dy, dt, c)
%       p_xy = kspaceLineRecon(p_ty, dy, dt, c, ...)
%
% INPUTS:
%       p_ty        - pressure time-series recorded over an evenly spaced
%                     array of sensor points on a line (indexed as t, y)
%       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.
%
%       'DataOrder' - String input which sets the data order (default =
%                     'ty'). Valid inputs are 'ty' and 'yt'
%       'Interp'    - string input controlling the interpolation method
%                     used by interp2 in the reconstruction (default =
%                     '*nearest')
%       'Plot'      - 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_xy        - estimate of the initial acoustic pressure
%                     distribution (indexed as x, y)
%
% ABOUT:
%       author      - Bradley Treeby and Ben Cox
%       date        - 11th January 2009
%       last update - 20th October 2011
%       
% This function is part of the k-Wave Toolbox (http://www.k-wave.org)
% Copyright (C) 2009-2012 Bradley Treeby and Ben Cox
%
% See also interp2, kspacePlaneRecon, makeGrid

% This file is part of k-Wave. k-Wave is free software: you can
% redistribute it and/or modify it under the terms of the GNU Lesser
% General Public License as published by the Free Software Foundation,
% either version 3 of the License, or (at your option) any later version.
% 
% k-Wave is distributed in the hope that it will be useful, but WITHOUT ANY
% WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
% FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for
% more details. 
% 
% You should have received a copy of the GNU Lesser General Public License
% along with k-Wave. If not, see <http://www.gnu.org/licenses/>.

% B.04 Modifications
% - removed dependence on makeGrid
% - modified inputs based on (z, x) to (x, y) re-ordering across k-Wave

% start timer
tic;

% define defaults
num_req_inputs = 4;
data_order = 'ty';
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 'DataOrder'
                data_order = varargin{input_index + 1};
                if ~strcmp(data_order, 'ty') && ~strcmp(data_order, 'yt')
                    error('Unknown setting for optional input DataOrder');
                end
            case 'Interp'
                interp_method = varargin{input_index + 1};
            case 'Plot'
                plot_recon = 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

% reorder the data if needed (p_ty)
if strcmp(data_order, 'yt')
    p = p.';
end

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

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

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

% create the computational grid for w; this represents the initial sampling
% of p(w, ky)
if rem(Nt, 2)
    % Nt is odd
    w = repmat(2*pi/dt*(-0.5:1/(Nt - 1):0.5).', 1, Ny);
else
    % Nt is even
    w = repmat(2*pi/dt*(-0.5:1/Nt:0.5 - 1/Nt).', 1, Ny);
end

% create the computational grid for ky
if rem(Ny, 2)
    % Nx is odd
    ky = repmat(2*pi/dy*(-0.5:1/(Ny - 1):0.5), Nt, 1);
else
    % Nx is even
    ky = repmat(2*pi/dy*(-0.5:1/Ny:0.5 - 1/Ny), Nt, 1);
end

% create the computational grid for kx
Nx = Nt;
dx = dt*c;
if rem(Nx, 2)
    % Nx is odd
    kx = repmat(2*pi/dx*(-0.5:1/(Nx - 1):0.5).', 1, Ny);
else
    % Nx is even
    kx = repmat(2*pi/dx*(-0.5:1/Nx:0.5 - 1/Nx).', 1, Ny);
end

% remap the computational grid for kx onto w using the dispersion
% relation w/c = sqrt(kx.^2 + ky.^2) 
kx = c.*sqrt( kx.^2 + ky.^2 );

% calculate the scaling factor using the value of kx, where
% kx = sqrt( (w/c).^2 - kgrid.ky.^2) and then manually replacing the DC
% value with its limit (otherwise NaN results)
sf = c^2*sqrt( (w/c).^2 - ky.^2 )./(2*w);
sf(w == 0 & ky == 0) = c/2;

% compute the FFT of the input data p(t, y) to yield p(w, ky) and scale
p = sf.*fftshift(fftn(fftshift(p)));

% remove unused variables
clear sf;

% exclude the inhomogeneous part of the wave
p(abs(w) < abs(c*ky)) = 0;

% compute the interpolation from p(w, ky) to p(kx, ky)and then force to be
% symmetrical 
p = interp2(ky, w, p, ky, kx, interp_method);

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

% compute the inverse FFT of p(kx, ky) to yield p(x, y)
p = real(ifftshift(ifftn(ifftshift(p))));

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

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

% enfore positivity condition
if positivity_cond
    disp('  applying positivity condition...');
    p(p < 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 (Nt/2)*dt*c]; 
    y_axis = [0 Ny*dy];
    
    % select suitable axis scaling factor
    [x_sc, scale, prefix] = scaleSI(max([Ny*dy, (Nt/2)*dt*c])); 
    
    % select suitable plot scaling factor
    plot_scale = max(p(:));
    
    % create the figure
    figure;
    imagesc(y_axis*scale, x_axis*scale, p, [-plot_scale, plot_scale]);
    axis image;
    colormap(getColorMap);
    xlabel(['Sensor Position [' prefix 'm]']);
    ylabel(['Depth [' prefix 'm]']);
    colorbar;
end