/* @(#)wa_deconv.c	17.1.1.1 (ES0-DMD) 01/25/02 17:21:23 */
/*===========================================================================
  Copyright (C) 1995 European Southern Observatory (ESO)
 
  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., 675 Massachusetss Ave, Cambridge, 
  MA 02139, USA.
 
  Corresponding concerning ESO-MIDAS should be addressed as follows:
	Internet e-mail: midas@eso.org
	Postal address: European Southern Observatory
			Data Management Division 
			Karl-Schwarzschild-Strasse 2
			D 85748 Garching bei Muenchen 
			GERMANY
===========================================================================*/

/****************************************************************************
**            Copyright (C) 1993 by European Southern Observatory
*****************************************************************************
**
**    UNIT
**
**    Version: 17.1.1.1
**
**    Author: Jean-Luc Starck
**
**    Date:  02/01/25
**    
**    File:  wa_deconv.c
**
*****************************************************************************
**
**    DESCRIPTION routines used for the deconvolution
**    -----------  
**
*****************************************************************************
**
** dec_pos_max (Tab, Nl, Nc, Ind_i, ind_j , Val_Max)
** float *Tab;
** int Nl, Nc;
** int *Ind_i, *ind_j;
** float *Val_Max;
**
** seek the maximum in an image
** Tab = image of size Nl*Nc
** Ind_i, ind_j = position of the max
** Val_Max = max values
**
*****************************************************************************
**
** dec_line_column (N, N_Out)
** int N, *N_Out;
**
** *N_Out =  the power of 2 superior or equal at N
**
*****************************************************************************
**
** int dec_test_ind (ind, N)
** int ind, N
**
** return the indice value with a miror at the bord
** if ind >= N val_return = 2 * (N - 1) - ind;
** if ind < 0 val_return = - ind
** if 0 <= ind < N val_return =  ind
**
*****************************************************************************
**
** dec_extract_ima (Imag, Nl1, Nc1, Imag_Out, Nl0, Nc0)
** float *Imag, *Imag_Out;
** int Nl0, Nc0, Nl1, Nc1;
**
** extract an image of size Nl0,Nc0 in a image Nl1,Nc1
** Imag_Out is the center of the image Imag
** we must have:
**            Nl0 < Nl1
**            Nc0 < Nc1
**
*****************************************************************************
**
** dec_insert_ima (Imag, Nl0, Nc0, Imag_Out, Nl1, Nc1)
** float *Imag, *Imag_Out;
** int Nl0, Nc0, Nl1, Nc1;
**
** insert an image Imag at the center of a  bigger one. 
** all pixels of Imag_Out are modified. A miror is used
** at the board
**
*****************************************************************************
**
** dec_center_psf (D_Beam, Nl_Beam, Nc_Beam, Psf, Nl, Nc)
** float *D_Beam, *Psf;
** int Nl, Nc, Nl_Beam, Nc_Beam;
**
** move the psf to the center of a bigger image
** 
*****************************************************************************
**
** dec_convol_conj (Resi, Psf_cf, Nl, Nc)
** complex_float *Psf_cf;
** float *Resi;
** int Nl,Nc;
**
** convolve the image Resi with the conjugate of the Fourier
** transform of the PSF
**
*****************************************************************************
** 
** dec_convol (Resi, Psf_cf, Result, Nl, Nc)
** complex_float *Psf_cf;
** float *Resi, *Result;
** int Nl,Nc;
**
** convolve the image Resi with the PSF
**
**
*****************************************************************************
** 
** dec_signif_struct(Resi,Nl0,Nc0,Nbr_Plan,Noise_Ima,N_Sigma,Type_Transform)
** float *Resi;
** int Nl0, Nc0;
** int Nbr_Plan, Type_Transform;
** float Noise_Ima, N_Sigma;
**
** suppress the noise in Resi
** Resi = image
** Nl0,Nc0 = image size
** Nbr_Plan = number of scales of the transform
** Noise_Ima = standard deviation of the noise in Resi
** N_Sigma = threshold parameter
** Type_Transform = chosen wavelet transform algorithm
**
*****************************************************************************
** 
** dec_wa_lucy (Imag, Obj, Resi, Psf_cf, Nl, Nc, Eps_cv, Noise_Ima, 
**              N_Sigma, Nbr_Plan, Nbr_Iter, Transform)
** float *Imag, *Obj, *Resi;
** complex_float *Psf_cf;
** int Nl,Nc,Nbr_Plan,Nbr_Iter, Transform;
** float Eps_cv, Noise_Ima, N_Sigma;
**
** deconvolution of an image by Lucy's algorithm with a regularization in 
** the wavelet space
**
** Imag = IN: data
** Obj = OUT: deconvolved object
** Resi = OUT: residual
** Psf_cf = IN: fourier transform of the PSF
** Nl,Nc = IN: image sizes
** Esp_cv = IN: convergence parameter
** Nbr_Plan = IN: number of scales of the transform
** Noise_Ima = IN: standard deviation of the noise in Resi
** N_Sigma = IN: threshold parameter
** Type_Transform = IN: chosen wavelet transform algorithm
** Nbr_Iter = IN: maximum number of iterations
** 
*****************************************************************************
** 
** dec_wa_make_psf (Psf, Nl, Nc, Fwhm)
** float *Psf;
** int Nl,Nc;
** float Fwhm;
**
** creates a gaussian with a Full Width at Half Maximum
** equal to Fwhm
**
*****************************************************************************
** 
** dec_wa_gradient (Imag, Obj, Resi, Psf_cf, Nl, Nc, Eps_cv, Noise_Ima, 
**                  N_Sigma, Nbr_Plan, Nbr_Iter, Transform)
** float *Imag, *Obj, *Resi;
** complex_float *Psf_cf;
** int Nl,Nc,Nbr_Plan,Nbr_Iter, Transform;
** float Eps_cv, Noise_Ima, N_Sigma;
** 
** deconvolution by the one-step gradient method
** with a regularization in the wavelet space
**
** Imag = IN: data
** Obj = OUT: deconvolved object
** Resi = OUT: residual
** Psf_cf = IN: fourier transform of the PSF
** Nl,Nc = IN: image sizes
** Esp_cv = IN: convergence parameter
** Nbr_Plan = IN: number of scales of the transform
** Noise_Ima = IN: standard deviation of the noise in Resi
** N_Sigma = IN: threshold parameter
** Type_Transform = IN: chosen wavelet transform algorithm
** Nbr_Iter = IN: maximum number of iterations
**
*****************************************************************************
** 
** dec_wa_cittert (Imag, Obj, Resi, Psf_cf, Nl, Nc, Eps_cv, Noise_Ima, 
**                N_Sigma, Nbr_Plan, Nbr_Iter, Transform, Fwmh)
** float *Imag, *Obj, *Resi;
** complex_float *Psf_cf;
** int Nl,Nc,Nbr_Plan,Nbr_Iter, Transform;
** float Eps_cv, Noise_Ima, N_Sigma, Fwmh;
**
** deconvolution by the van citter method
** with a regularization in the wavelet space
** if Fwmh > 0 then the resolution is limited. 
**
** Imag = IN: data
** Obj = OUT: deconvolved object
** Resi = OUT: residual
** Psf_cf = IN: fourier transform of the PSF
** Nl,Nc = IN: image sizes
** Esp_cv = IN: convergence parameter
** Nbr_Plan = IN: number of scales of the transform
** Noise_Ima = IN: standard deviation of the noise in Resi
** N_Sigma = IN: threshold parameter
** Type_Transform = IN: chosen wavelet transform algorithm
** Nbr_Iter = IN: maximum number of iterations
** Fwmh = IN: limited resolution paramter
**
*****************************************************************************/

static char sccsid[] = "@(#)wa_deconv.c 17.1.1.1 02/01/25 ESO 1993 @(#)";

#include <stdio.h>
#include <math.h>
#include <string.h>

#include<midas_def.h>

#include "Def_Math.h"
#include "Def_Mem.h"
#include "Def_Wavelet.h"

#define VISU_PARAM TRUE
#define VISU_DATA TRUE

/****************************************************************************/

dec_pos_max (Tab, Nl, Nc, Ind_i, ind_j , Val_Max)
float *Tab;
int Nl, Nc;
int *Ind_i, *ind_j;
float *Val_Max;
{
   int i,j;
   float Val_Abs_Max;
   int Bande = 1;

   *Val_Max = 0.;
   Val_Abs_Max = 0.;
   for (i = Bande; i < Nl-Bande; i++)
   {
      for (j = Bande; j < Nc-Bande; j++)
      {
         if (Tab [i*Nc+j] > Val_Abs_Max)
         {
            Val_Abs_Max = Tab [i*Nc+j];
            *Val_Max = Tab [i*Nc+j];
            *Ind_i = i;
            *ind_j = j;
         }
      }
   }
}

/***************************************************************************/

dec_line_column (N, N_Out)
int N, *N_Out;
/* We must check Nl to make sure it is a power of 2 */
{
    int len_exp,temp;

    len_exp = (int)(0.3+log((double)(N))/(log(2.0)));
    INT_POW(2,len_exp,temp);
    if (temp < N) temp *= 2; 
    *N_Out = temp;
}

/**************************************************************************/

int dec_test_ind (ind, N)
int ind, N;
{
    int Val;

    if (ind < 0) Val = - ind;
    else
    {
        if (ind >= N) Val = 2 * (N - 1) - ind;
        else Val = ind;
    }
    if ((Val >= N) || (Val < 0)) Val = -1;
    return (Val);
}

/**************************************************************************/

dec_extract_ima (Imag, Nl1, Nc1, Imag_Out, Nl0, Nc0)
float *Imag, *Imag_Out;
int Nl0, Nc0, Nl1, Nc1;
{
    int i1,j1,i0,j0, Depi, Depj;

    Depi = (Nl1 - Nl0) / 2;
    Depj = (Nc1 - Nc0) / 2;

    for (i0 = 0; i0 < Nl0; i0++)
    for (j0 = 0; j0 < Nc0; j0++)
    {
       i1 = i0 + Depi;
       j1 = j0 + Depj;
       Imag_Out [i0*Nc0+j0] = Imag [i1*Nc1+j1];
    }
}

/**************************************************************************/

dec_insert_ima (Imag, Nl0, Nc0, Imag_Out, Nl1, Nc1)
float *Imag, *Imag_Out;
int Nl0, Nc0, Nl1, Nc1;
{
    int i1,j1,i0,j0, Depi, Depj;

    Depi = (Nl1 - Nl0) / 2;
    Depj = (Nc1 - Nc0) / 2;

    for (i1 = 0; i1 < Nl1; i1++)
    for (j1 = 0; j1 < Nc1; j1++)
    {
       i0 = dec_test_ind (i1 - Depi, Nl0);
       j0 = dec_test_ind (j1 - Depj, Nc0);
       if ((i0 < 0) || (j0 < 0)) Imag_Out [i1*Nc1+j1] = 0.;
       else Imag_Out [i1*Nc1+j1] = Imag [i0*Nc0+j0];
    }
}

/*********************************************************************/

dec_center_psf (D_Beam, Nl_Beam, Nc_Beam, Psf, Nl, Nc)
float *D_Beam, *Psf;
int Nl, Nc, Nl_Beam, Nc_Beam;
{
   int Ind_i, Ind_j,i,j,Dep_i, Dep_j;
   float Val_Max;

   dec_pos_max (D_Beam, Nl_Beam, Nc_Beam, &Ind_i, &Ind_j , &Val_Max);

   for (i=0;i<Nl*Nc;i++) Psf[i] = 0.;

   for (i = 0; i < Nl_Beam; i++)
   for (j = 0; j < Nc_Beam; j++)
   {
      Dep_i = Ind_i - i + Nl / 2;
      Dep_j = Ind_j - j + Nc / 2;
      if ((Dep_i >= 0) && (Dep_i < Nl) && (Dep_j >= 0) && (Dep_j < Nc))
      {
         Psf[Dep_i*Nc+Dep_j] = D_Beam[i*Nl_Beam+j];
      }
    }
}

/****************************************************************************/

dec_convol_conj (Resi, Psf_cf, Nl, Nc)
complex_float *Psf_cf;
float *Resi;
int Nl,Nc;
{
   complex_float *O_cf, Val,Val_Psf;
   int j;

   O_cf = cf_vector_alloc (Nl*Nc);
   prepare_fft_real (Resi, O_cf, Nl);
   ft_cf_any_power_of_2 (O_cf, 1,  Nl);

   for (j = 0; j < Nl*Nc; j++) 
   {
      Val_Psf.re = Psf_cf[j].re;
      Val_Psf.im = - Psf_cf[j].im;
      CF_MLT(Val_Psf, O_cf[j], Val);
      O_cf[j].re = Val.re;
      O_cf[j].im = Val.im;
   }

   ft_cf_any_power_of_2 (O_cf, -1,  Nl);
   for (j = 0; j < Nl*Nc; j++) Resi[j] = O_cf[j].re;

   free ((char *) O_cf);
}

/***************************************************************************/

dec_convol (Resi, Psf_cf, Result, Nl, Nc)
complex_float *Psf_cf;
float *Resi, *Result;
int Nl,Nc;
{
   complex_float *O_cf, Val;
   int j;

   O_cf = cf_vector_alloc (Nl*Nc);
   prepare_fft_real (Resi, O_cf, Nl);
   ft_cf_any_power_of_2 (O_cf, 1,  Nl);

   for (j = 0; j < Nl*Nc; j++) 
   {
       CF_MLT(Psf_cf[j], O_cf[j], Val);
       O_cf[j].re = Val.re;
       O_cf[j].im = Val.im;
   }

   ft_cf_any_power_of_2 (O_cf, -1,  Nl);
   for (j = 0; j < Nl*Nc; j++) Result[j] = O_cf[j].re;

   free ((char *) O_cf);
}

/************************************************************************/

dec_signif_struct (Resi, Nl0, Nc0, Nbr_Plan, Noise_Ima,N_Sigma,Type_Transform)
float *Resi;
int Nl0, Nc0;
int Nbr_Plan, Type_Transform;
float Noise_Ima, N_Sigma;
{
   int i,j,Nl,Nc;
   wave_transf_des W_Resi;
   float *Plan_Resi;
   float Noise;
   float Coef,Fc = 0.5, lib_mat_ecart_type ();

   wavelet_transform_data (Resi,Nl0, Nc0,&W_Resi,Type_Transform, Fc, Nbr_Plan);

   for (i = 0; i < Nbr_Plan-1; i++)
   {
       wavelet_pointer_plan (&W_Resi, &Plan_Resi, &Nl, &Nc, i+1, 0);
            
       switch(i)
       {
          case 0: Noise = 0.89*Noise_Ima;Coef=N_Sigma;break;
          case 1: Noise = 0.2*Noise_Ima;Coef=N_Sigma;break;
	  case 2: Noise = 0.086*Noise_Ima;Coef=N_Sigma;break;
	  case 3: Noise = 0.04*Noise_Ima;Coef=N_Sigma;break;
          default: Noise = 0.;
       }
       Noise *= Coef;

       for (j = 0; j < Nl*Nc; j++)
           if (fabs(Plan_Resi[j]) < Noise)  Plan_Resi[j] = 0.;
   }
   wavelet_reconstruct_data (&W_Resi, Resi, 1);
   wave_io_free (&W_Resi);
}

/*********************************************************************/

dec_wa_lucy (Imag, Obj, Resi, Psf_cf, Nl, Nc, Eps_cv, Noise_Ima, 
              N_Sigma, Nbr_Plan, Nbr_Iter, Transform)
float *Imag, *Obj, *Resi;
complex_float *Psf_cf;
int Nl,Nc,Nbr_Plan,Nbr_Iter, Transform;
float Eps_cv, Noise_Ima, N_Sigma;
{
    int j,Iter=0, Size = Nl*Nc;
    float lib_mat_ecart_type();
    float Moy, Sigma, Old_Sigma, Delta;
    float *Imag_n, Noise;
    char Send[200];

    Imag_n = f_vector_alloc (Nl*Nc);

    /* Noise estimation 
       if Noise_Ima = 0 then no values has been introduced
       by the user and the program has to estimate the noise 
    */
    if (Noise_Ima < FLOAT_EPSILON)
            lib_mat_detect_snr (Nc,Nl,Imag,1,3,&Moy,&Noise);
    else Noise = Noise_Ima;
 
    /* We choose the filtered image as the first
       estimation of the object */
    wave_filter_imag (Imag, Nl, Nc, Obj, 5., FILTER_TRESHOLD,
                      1, 2, Nbr_Plan, 0.5, Noise);

    /* The first estimation has to be > 0 */
    for (j = 0; j < Size; j++) if (Obj[j] < 0.) Obj[j] = 0.;
 
    /* Initialization */
    Delta = Sigma = 1.e20;
    do
    {
       Old_Sigma = Sigma; 

       /* Residual estimation */ 
       dec_convol (Obj, Psf_cf, Imag_n, Nl, Nc);
       for (j = 0; j < Size; j++) Resi[j] = Imag[j] - Imag_n[j];

       /* Standard deviation of the residual */
       lib_mat_moy_ecart_type (Resi, Nl, Nc, &Sigma, &Moy);

       /* New estimation of the noise */
       if ((Noise > Sigma) && (Delta > 0.01)) Noise = Sigma;

       /* Significant structure extraction */
       dec_signif_struct (Resi, Nl, Nc, Nbr_Plan, Noise, N_Sigma, Transform);

       /* Calculate the multiplicate term in Lucy's iteration */
       for (j = 0; j < Size; j++) 
       {
          if (fabs(Imag_n[j]) > FLOAT_EPSILON)
                 Resi[j] = (Imag_n[j] + Resi[j]) / Imag_n[j];
          else Resi[j] = 1.;
          if (Resi[j] < 0.)  Resi[j] = 1.;

       }
       dec_convol_conj (Resi, Psf_cf, Nl, Nc); 

       /* calculate the next estimation of the object */
       for (j = 0; j < Size; j++) Obj[j] *= Resi[j];

       Delta = (Old_Sigma-Sigma)/Sigma;

#if VISU_DATA
       if ((Iter > 0) && (Iter % 1 == 0))
       {
          sprintf(Send,"%d: Sigma, Average residual : %f, %f",Iter,Sigma,Moy);
          SCTPUT(Send);
          sprintf (Send, "   Cvg parameter: %f", Delta);
          SCTPUT(Send);
       }
#endif
       Iter ++;
    } while ((Iter < Nbr_Iter) && (Delta > Eps_cv));

   /* Residual */
   dec_convol (Obj, Psf_cf, Imag_n, Nl, Nc);
   for (j = 0; j < Size; j++) Resi[j] = Imag[j] - Imag_n[j];

   free ((char *) Imag_n);
}

/****************************************************************************/

dec_wa_make_psf (Psf, Nl, Nc, Fwhm)
float *Psf;
int Nl,Nc;
float Fwhm;
{
    int Demi_Size,Dist,Delta_i;
    float sigma,sigma2,Energ = 0.;
    register int i,j;

    Demi_Size = Nl / 2;

    sigma = 0.5 * Fwhm / sqrt (2. * log ((double) 2.));
    sigma2 = -2. * sigma * sigma;

    for (i = 0; i < Nl; i++)
    {
        Delta_i = (i - Demi_Size) * (i - Demi_Size);
        for (j = 0; j < Nc; j++)
	{
	    Dist = Delta_i + (j - Demi_Size) * (j - Demi_Size);
	    Psf [i*Nc+j] = exp((double)((float) Dist / sigma2));
	    Energ +=  Psf [i*Nc+j];
	 }
    } 

    for (i = 0; i < Nl*Nc; i++) Psf [i] /= Energ;
}

/***************************************************************************/

dec_wa_gradient (Imag, Obj, Resi, Psf_cf, Nl, Nc, Eps_cv, Noise_Ima, 
                 N_Sigma, Nbr_Plan, Nbr_Iter, Transform)
float *Imag, *Obj, *Resi;
complex_float *Psf_cf;
int Nl,Nc,Nbr_Plan,Nbr_Iter, Transform;
float Eps_cv, Noise_Ima, N_Sigma;
{
    int j,Iter=0, Size = Nl*Nc;
    float lib_mat_ecart_type();
    float Moy, Sigma, Old_Sigma, Delta;
    float *Imag_n, Noise;
    char Send[200];

    Imag_n = f_vector_alloc (Nl*Nc);

    /* Noise estimation 
       if Noise_Ima = 0 then no values has been introduced
       by the user and the program has to estimate the noise 
    */
    if (Noise_Ima < FLOAT_EPSILON)
            lib_mat_detect_snr (Nc,Nl,Imag,1,3,&Moy,&Noise);
    else Noise = Noise_Ima;
 
    /* We choose the filtered image as the first
       estimation of the object */
    wave_filter_imag (Imag, Nl, Nc, Obj, 5., FILTER_TRESHOLD,
                      1, 2, Nbr_Plan, 0.5, Noise);

    /* The first estimation has to be > 0 */
    for (j = 0; j < Size; j++) if (Obj[j] < 0.) Obj[j] = 0.;
 
    /* Initialization */
    Delta = Sigma = 1.e20;
    do
    {
       Old_Sigma = Sigma; 

       /* Residual estimation */ 
       dec_convol (Obj, Psf_cf, Imag_n, Nl, Nc);
       for (j = 0; j < Size; j++) Resi[j] = Imag[j] - Imag_n[j];

       /* Standard deviation of the residual */
       lib_mat_moy_ecart_type (Resi, Nl, Nc, &Sigma, &Moy);

       /* New estimation of the noise */
       if ((Noise > Sigma) && (Delta > 0.01)) Noise = Sigma;

       /* Significant structure extraction */
       dec_signif_struct (Resi, Nl, Nc, Nbr_Plan, Noise, N_Sigma, Transform);

       /* Convolves with conjugate of the PSF */
       dec_convol_conj (Resi, Psf_cf, Nl, Nc); 

       /* calculate the next estimation of the object */
       for (j = 0; j < Size; j++)
       { 
           Obj[j] += Resi[j];

           /* Positivity */
           if (Obj[j] < 0.) Obj[j] = 0.;
       }

       Delta = (Old_Sigma-Sigma)/Sigma;

#if VISU_DATA
       if ((Iter > 0) && (Iter % 1 == 0))
       {
          sprintf(Send,"%d: Sigma, Average residual : %f, %f",Iter,Sigma,Moy);
          SCTPUT(Send);
          sprintf (Send, "   Cvg parameter: %f", Delta);
          SCTPUT(Send);
       }
#endif
       Iter ++;
    } while ((Iter < Nbr_Iter) && (Delta > Eps_cv));

   /* Residual */
   dec_convol (Obj, Psf_cf, Imag_n, Nl, Nc);
   for (j = 0; j < Size; j++) Resi[j] = Imag[j] - Imag_n[j];

   free ((char *) Imag_n);
}

/****************************************************************************/

static cittert_mult_cf (Imag1, Imag2, Imag3, Nl, Nc)
complex_float *Imag1, *Imag2, *Imag3;
int Nl, Nc;
{
   int j;
   complex_float Val;

   for (j = 0; j < Nl*Nc; j++) 
   {
      CF_MLT(Imag1[j], Imag2[j], Val);
      CF_ASS(Val, Imag3[j]);
   }
}

/****************************************************************************/

dec_wa_cittert (Imag, Obj, Resi, Psf_cf, Nl, Nc, Eps_cv, Noise_Ima, 
                N_Sigma, Nbr_Plan, Nbr_Iter, Transform, Fwmh)
float *Imag, *Obj, *Resi;
complex_float *Psf_cf;
int Nl,Nc,Nbr_Plan,Nbr_Iter, Transform;
float Eps_cv, Noise_Ima, N_Sigma, Fwmh;
{
    int j,Iter=0, Size = Nl*Nc;
    float lib_mat_ecart_type();
    float Moy, Sigma, Old_Sigma, Delta;
    float *Imag_n, Noise;
    char Send[200];
    float *Regul;
    complex_float *Regul_cf;

    Imag_n = f_vector_alloc (Nl*Nc);
    
    if (Fwmh > FLOAT_EPSILON)
    {
        Regul  = f_vector_alloc (Nl*Nc);
        Regul_cf = cf_vector_alloc (Nl*Nc);

        dec_wa_make_psf (Regul, Nl, Nc, Fwmh);
        prepare_fft_real (Regul, Regul_cf, Nl);
        ft_cf_any_power_of_2 (Regul_cf, 1,  Nl);
        free ((char *) Regul);

        /* convolution of the PSF by the regularization term */
        cittert_mult_cf (Psf_cf, Regul_cf, Psf_cf, Nl, Nc);
    }

    /* Noise estimation 
       if Noise_Ima = 0 then no values has been introduced
       by the user and the program has to estimate the noise 
    */
    if (Noise_Ima < FLOAT_EPSILON)
            lib_mat_detect_snr (Nc,Nl,Imag,1,3,&Moy,&Noise);
    else Noise = Noise_Ima;
 
    /* We choose the filtered image as the first
       estimation of the object */
    wave_filter_imag (Imag, Nl, Nc, Obj, 5., FILTER_TRESHOLD,
                      1, 2, Nbr_Plan, 0.5, Noise);

    /* The first estimation has to be > 0 */
    for (j = 0; j < Size; j++) if (Obj[j] < 0.) Obj[j] = 0.;
 
    /* Initialization */
    Delta = Sigma = 1.e20;
    do
    {
       Old_Sigma = Sigma; 

       /* Residual estimation */ 
       dec_convol (Obj, Psf_cf, Imag_n, Nl, Nc);
       for (j = 0; j < Size; j++) Resi[j] = Imag[j] - Imag_n[j];

       /* Standard deviation of the residual */
       lib_mat_moy_ecart_type (Resi, Nl, Nc, &Sigma, &Moy);

       /* New estimation of the noise */
       if ((Noise > Sigma) && (Delta > 0.01)) Noise = Sigma;

       /* Significant structure extraction */
       dec_signif_struct (Resi, Nl, Nc, Nbr_Plan, Noise, N_Sigma, Transform);

       /* calculate the next estimation of the object */
       for (j = 0; j < Size; j++)
       {
           Obj[j] += Resi[j];

           /* Positivity */
           if (Obj[j] < 0.) Obj[j] = 0.;
       }
       Delta = (Old_Sigma-Sigma)/Sigma;

#if VISU_DATA
       if ((Iter > 0) && (Iter % 1 == 0))
       {
          sprintf(Send,"%d: Sigma, Average residual : %f, %f",Iter,Sigma,Moy);
          SCTPUT(Send);
          sprintf (Send, "   Cvg parameter: %f", Delta);
          SCTPUT(Send);
       }
#endif
       Iter ++;
    } while ((Iter < Nbr_Iter) && (Delta > Eps_cv));

   /* Residual */
   dec_convol (Obj, Psf_cf, Imag_n, Nl, Nc);
   for (j = 0; j < Size; j++) Resi[j] = Imag[j] - Imag_n[j];

   if (Fwmh > FLOAT_EPSILON)
   {
       /* limited resolution object */
       dec_convol (Obj, Regul_cf, Obj, Nl, Nc);
       free ((char *) Regul_cf);
   }
   free ((char *) Imag_n);
}

/****************************************************************************/