function psf = LsfToPsf(lsf,varargin)
%LSFToPSF  Convert a line spread function to a point spread function.
%   psf = LSFTOPSF(lsf,varargin)
%
%   This works by taking the lsf into the one-dimensional frequency domain
%   (i.e. get the 1D MTF) and then creating a circularly symmetric version
%   of this.  The 2D frequency representation is then converted back to the
%   spatial domain to produce the psf.  This produces a spatially-symmetric
%   PSF consistent with the measured line spread function.
%
%   This method is described by Marchand, 1964, JOSA, 54, 7, pp. 915-919
%   and is one of several methods provided.  In 1964, taking the Fourier
%   transform was computationally intense. 
%
%   There is a second paper by Marchand (1965, JOSA, 55, 4, 352-354) which
%   treats the more general case where you have line spread functions for
%   many orientations and want to recover an psf that is not necessarily
%   spatially-symmetric.
%
%   The lsf must be spatially symmetric.  This makes sense given that we
%   are going to recover a spatially symmetric psf.
%
%   You want to make sure that the spatial suport is large enough to
%   capture the full lsf.
%
%   The returned psf is normalized to have unit volume.
%
%   See also PSFTOLSF, PSYCHOPTICSTEST

% Handle even versus odd dimension of spatial support.
n1D = length(lsf);
if (mod(n1D,2) ~= 0)
    centerPosition = floor(n1D/2) + 1;
    negValues = lsf(centerPosition-1:-1:1);
    posValues = lsf(centerPosition+1:end);
else
    centerPosition = n1D/2 + 1;
    negValues = lsf(centerPosition-1:-1:2);
    posValues = lsf(centerPosition+1:end);
end
if (any(posValues ~= negValues))
    error('LSF is not spatially symmetric');
end
    
% Convert LSF to 1-D MTF
%
% Taking the absolute value.  The fftshift calls wrapping the fft lets us
% work in a coordinate system where 0 in both space and frequency is at
% position n1D/2+1.  If n1D were odd, we'd have to figure out more
% precisely the difference between fftshift and ifftshift, this may not be
% used completely correctly here.  My current view is that ifftshift takes
% one from 0 centered to 0 at first entry, and fftshift goes the other way.
lsfNormalized = lsf/sum(lsf(:));
lsfMTFCenteredRaw = fftshift(fft(ifftshift(lsf)));

% Result should be real.  Check that it's OK to tolerance and then get rid
% of any pesky imaginary components.
lsfOTFCenteredImag = imag(lsfMTFCenteredRaw);
if (any(abs(lsfOTFCenteredImag > 1e-10)))
    error('Imaginary part of 1D MTF is bigger than makes sense from numerical error alone');
end
lsfOTFCentered = abs(lsfMTFCenteredRaw);

% Create symmetric 2D MTF from 1D MTF
%
% Values outside the frequency support of the lsf-based OTF are set to
% zero.  This should be fine as long as the spatial support is large enough
% so that the MTF has fallen to zero at highest frequencies specified.
%
% Need to handle even versus odd spatial support here.
radiusMatrix = MakeRadiusMat(n1D,n1D,centerPosition);
if (rem(n1D,2) == 0)
    psfOTFCentered = interp1(0:n1D/2-1,lsfOTFCentered(centerPosition:end),radiusMatrix,'linear',0);
else
    psfOTFCentered = interp1(0:floor(n1D/2),lsfOTFCentered(centerPosition:end),radiusMatrix,'linear',0);
end

% Do inverse fft, now in 2D.  Again wrap with ifftshift and fftshift so that we
% are working in 0 centered coordinates.
psfRaw = fftshift(ifft2(ifftshift(psfOTFCentered)));

% Result should be real.  Check that it's OK to tolerance and then get rid
% of any pesky imaginary components.
psfImag = imag(psfRaw);
if (any(abs(psfImag) > 1e-10))
    error('Imaginary component of derived psf too big for comfort.');
end
psf = abs(psfRaw);

% Normalize to unit volume
psf = psf/sum(psf(:));

end