/*
   LUCY.I
   Lucy-Richardson-type image deconvolution

    Copyright (C) 2008  Thibaut Paumard <paumard@users.sourceforge.net>

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program 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 General Public License for more details.

    You should have received a copy of the GNU General Public License along
    with this program; if not, write to the Free Software Foundation, Inc.,
    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

   History:

   v. 0.1, Fri, 08 Feb 2008 09:17:04 +0100 Thibaut Paumard
      initial release.

*/

extern lucy_convolve;
/* DOCUMENT convolved=lucy_convolve(image,psf)

     convolve an IMAGE by a point spread function (PSF).

     lucy_convolve must be defined for lucy() to work, but the lucy
     package does not provide it. By default, If yeti_fftw is
     installed, lucy_convolve is a copy of fftw_convolve. Else, if
     yutils is installed, lucy_convolve is the same as
     fft_convolve. If neither of these two packages is installed, it's
     up to you to provide a compatible convolution routine.

     If you install either yutils or yeti_fftw while Yorick is
     running, please re-include lucy.i so that the right FFT library
     is detected.

   SEE ALSO: lucy, fftw_convolve, fft_convolve
 */

if (!is_func(lucy_convolve)) {
  if (numberof(find_in_path("yeti_fftw.i"))) {
    require, "yeti_fftw.i";
    lucy_convolve=fftw_convolve;
  } else if (numberof(find_in_path("fft_utils.i"))) {
    require, "fft_utils.i";
    lucy_convolve=fft_convolve;
  }    
}

func lucy(input,psf,&output,maxiter=) {
/* DOCUMENT deconv=lucy(image,psf [,guess])

     Deconvolve IMAGE using the Lucy - Richardson algorithm.

     Image and PSF must have the same shape. An optional GUESS can be
     provided (for instance the output from a previous run, to iterate
     some more). Keyword MAXITER controls the number of iterations
     (defaults to 100).

   SEE ALSO: clean
 */
  if (!is_func(lucy_convolve))
    error,"lucy_convolve() is not defined. Read help,lucy_convolve."
  if (is_void(maxiter)) maxiter=100;
  if (is_void(output)) output=input;
  for (n=1; n<=maxiter; n++) {
    reconv=lucy_convolve(output,psf);
    truc=input/reconv;
    output=output*lucy_convolve(truc,psf(0:1:-1,0:1:-1));
  }
  return output;
}