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