% Richardson-Lucy deconvolution (1D)

function result = rl(data,psf,Niter)

% Input arguments:
% data  = row vector of data to be deconvolved
% psf   = row vector of known point spread function (centered on element 1).
%         [Note: Internally, rl imposes unitarity on the FFT of psf.]
% Niter = number of Richardson-Lucy iterations to be performed
%
% CCK 2016-Apr-25


% Some preparatory calculatons...
psf_f = fft(psf);
psf_f = psf_f ./ psf(1); % normalize
flip_psf = fliplr(psf);
flip_psf = shift(flip_psf, 1); 
   % This line fixes the flipped PSF so that the central element 
   % (which lives at index 1) stays there upon flipping. If I
   % omit this shift, then RL fails to converge!
flip_psf_f = fft(flip_psf);
flip_psf_f = flip_psf_f ./ psf(1); % normalize

% Richardson-Lucy iterations...
result = data; % initial guess
for t=1:Niter
   correction_f = fft( data ./ ifft(fft(result) .* psf_f) ) .* flip_psf_f;
   correction = real( ifft(correction_f) );
   result = result .* correction;
end