/* @(#)wa_grad.c 17.1.1.1 (ES0-DMD) 01/25/02 17:21:40 */ /*=========================================================================== 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_grad.c ** ******************************************************************************* ** ** DESCRIPTION deconvolution by the one-step gradient method ** ----------- with a regularization in the wavelet space. ** ** ** PARAMETRES ** ---------- ** Imag: (keyword: IN_A): input image ** ** PSF: (keyword: IN_B): input PSF ** ** Object: (keyword: OUT_A): output deconvolved object ** ** Residual: (keyword: OUT_B): output residual ** ** Nbr_Plan: (keyword: INPUTI[1]): number of scales ** it is not necessary to take values for Nbr_Plan superior to 5 ** ** N_Sigma: (keyword: INPUTR[1]): If we threshold, the level is estimated b: ** Level = N_Sigma * Noise_Standart_Deviation ** N_Sigma = 3 is a standart value ** ** Sigma_Noise: (keyword: INPUTR[2]): ** standard deviation of the noise ** if Sigma_Noise = 0 then the standard deviation of the noise ** ** Eps_cv: (keyword: INPUTR[3]): convergence parameter ** ** Nbr_Iter: (keyword: INPUTI[2]): ** maximum number of iterations ** ** RESULTS ** ------- ** two files are created. The first contains the deconvolved object ** and the second the residual ** ******************************************************************************/ static char sccsid[] = "@(#)wa_grad.c 17.1.1.1 02/01/25 ESO 1993 @(#)"; #include #include #include #include #include "Def_Math.h" #include "Def_Mem.h" #include "Def_Wavelet.h" #define VISU_PARAM TRUE /****************************************************************************/ main() { int Nl,Nc,Nl_Psf,Nc_Psf,Nl1,Nc1; int Nbr_Plan, Nbr_Iter; string Name_File_Imag, Name_File_Psf, Name_File_Out, Name_File_Resi; float *Resi, *Psf, *Imag, *Imag1, *Psf1, *Result; float Eps_cv, N_Sigma, Noise_Ima; complex_float *Psf_cf; int Unit; int Null, Actvals, Buffer_Int, Maxvals, Felem; int Stat, Dummy; int Type_Transform = 2; /* Initialisation */ SCSPRO("wa_grad"); /* read the image 1 file name */ Felem = 1; Maxvals = 60; Stat = SCKGETC("IN_A", Felem, Maxvals, &Actvals, Name_File_Imag); /* read the PSF file name */ Felem = 1; Maxvals = 60; Stat = SCKGETC("IN_B", Felem, Maxvals, &Actvals, Name_File_Psf); /* read the number of scales */ Felem = 1; Maxvals = 1; Stat = SCKRDI("INPUTI", Felem, Maxvals, &Actvals, &Buffer_Int, &Unit, &Null); Nbr_Plan = Buffer_Int; /* read the maximum iteration number */ Felem = 2; Maxvals = 1; Stat = SCKRDI("INPUTI", Felem, Maxvals, &Actvals, &Buffer_Int, &Unit, &Null); Nbr_Iter = Buffer_Int; /* read the N_Sigma_Parameter */ Felem = 1; Maxvals = 1; Stat = SCKRDR("INPUTR", Felem, Maxvals, &Actvals, &N_Sigma, &Unit, &Null); /* read the noise standard deviation */ Felem = 2; Stat = SCKRDR("INPUTR", Felem, Maxvals, &Actvals, &Noise_Ima,&Unit, &Null); /* read the Convergence Parameter */ Felem = 3; Stat = SCKRDR("INPUTR", Felem, Maxvals, &Actvals, &Eps_cv, &Unit, &Null); /* read the output file name */ Felem = 1; Maxvals = 60; Stat = SCKGETC("OUT_A", Felem, Maxvals, &Actvals, Name_File_Out); /* read the residual file name */ Felem = 1; Stat = SCKGETC("OUT_B", Felem, Maxvals, &Actvals, Name_File_Resi); /* read the input images */ io_read_file_to_pict_f (Name_File_Imag, &Imag, &Nl, &Nc); io_read_file_to_pict_f (Name_File_Psf, &Psf, &Nl_Psf, &Nc_Psf); #if VISU_PARAM { char Send[100]; sprintf(Send," File_Name_Imag = %s\n", Name_File_Imag); SCTPUT(Send); sprintf(Send," File_Name_Psf = %s\n", Name_File_Psf); SCTPUT(Send); sprintf(Send," File_Name_Out = %s\n", Name_File_Out); SCTPUT(Send); sprintf(Send," File_Name_Resi = %s\n", Name_File_Resi); SCTPUT(Send); sprintf(Send," Nbr_Plan = %d\n", Nbr_Plan); SCTPUT(Send); sprintf(Send," Max Iter = %d\n", Nbr_Iter); SCTPUT(Send); sprintf(Send," Noise = %f\n", Noise_Ima); SCTPUT(Send); sprintf(Send," N_Sigma = %f\n", N_Sigma); SCTPUT(Send); sprintf(Send," Eps = %f\n", Eps_cv); SCTPUT(Send); sprintf(Send," File_Name_Out = %s\n", Name_File_Out); SCTPUT(Send); sprintf(Send," File_Name_Resi = %s\n", Name_File_Resi); SCTPUT(Send); } #endif /* the image is inserted in a bigger square image with a size which is a power of two */ dec_line_column (Nl, &Nl1); dec_line_column (Nc, &Nc1); if (Nl1 < Nc1) Nl1 = Nc1; else Nc1 = Nl1; #if VISU_PARAM printf ("Image Size: %d\n", Nl1); #endif Imag1 = f_vector_alloc (Nl1*Nc1); dec_insert_ima (Imag, Nl, Nc, Imag1, Nl1, Nc1); /* Compute the FFT of the PSF in an image with the same size */ Psf1 = f_vector_alloc (Nl1*Nc1); dec_center_psf (Psf, Nl_Psf, Nc_Psf, Psf1, Nl1, Nc1); NORM_ENERG(Psf1, Nl1*Nc1); Psf_cf = cf_vector_alloc (Nl1*Nc1); prepare_fft_real (Psf1, Psf_cf, Nl1); ft_cf_any_power_of_2 (Psf_cf, 1, Nl1); free ((char *) Psf); free ((char *) Psf1); Resi = f_vector_alloc (Nl1*Nc1); Result = f_vector_alloc (Nl1*Nc1); dec_wa_gradient (Imag1, Result, Resi, Psf_cf, Nl1, Nc1, Eps_cv, Noise_Ima, N_Sigma, Nbr_Plan, Nbr_Iter, Type_Transform); dec_extract_ima (Result, Nl1, Nc1, Imag, Nl, Nc); io_write_pict_f_to_file (Name_File_Out, Imag, Nl, Nc); dec_extract_ima (Resi, Nl1, Nc1, Imag, Nl, Nc); io_write_pict_f_to_file (Name_File_Resi, Imag, Nl, Nc); free ((char *) Resi); free ((char *) Imag); free ((char *) Imag1); free ((char *) Psf_cf); free ((char *) Result); SCSEPI(); } /***************************************************************************/