From f717c221e532655b184f30bc374fdc9ec251a1c5 Mon Sep 17 00:00:00 2001 From: Samuel Farrens Date: Thu, 4 Jul 2019 17:08:14 +0200 Subject: [PATCH 1/7] added mr1d files --- CMakeLists.txt | 4 +- src/im1d_convert.cc | 123 +++ src/im1d_der_loggamma.cc | 608 +++++++++++++++ src/im1d_info.cc | 201 +++++ src/im1d_period.cc | 398 ++++++++++ src/im1d_stf.cc | 258 ++++++ src/im1d_tend.cc | 152 ++++ src/im1d_tfreq.cc | 252 ++++++ src/libsparse1d/MR1D_Deconv.cc | 1336 ++++++++++++++++++++++++++++++++ src/libsparse1d/MR1D_Deconv.h | 169 ++++ src/mr1d_deconv.cc | 564 ++++++++++++++ src/mr1d_detect.cc | 332 ++++++++ src/mr1d_filter.cc | 679 ++++++++++++++++ src/mr1d_trans.cc | 414 ++++++++++ 14 files changed, 5489 insertions(+), 1 deletion(-) create mode 100644 src/im1d_convert.cc create mode 100644 src/im1d_der_loggamma.cc create mode 100644 src/im1d_info.cc create mode 100644 src/im1d_period.cc create mode 100644 src/im1d_stf.cc create mode 100644 src/im1d_tend.cc create mode 100644 src/im1d_tfreq.cc create mode 100755 src/libsparse1d/MR1D_Deconv.cc create mode 100755 src/libsparse1d/MR1D_Deconv.h create mode 100644 src/mr1d_deconv.cc create mode 100644 src/mr1d_detect.cc create mode 100644 src/mr1d_filter.cc create mode 100644 src/mr1d_trans.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index c864c89..473880e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -144,7 +144,9 @@ ${FFTW_FLAGS}" # Compile and link executables set(BINMR2D mr_transform mr_recons mr_filter mr_deconv mr_info cur_contrast - cur_deconv cur_filter cur_stat cur_trans) + cur_deconv cur_filter cur_stat cur_trans im1d_convert im1d_info im1d_stf + im1d_tfreq mr1d_detect mr1d_trans im1d_der_loggamma im1d_period im1d_tend + mr1d_filter) if(SPARSE3D) set(BINMR2D ${BINMR2D} mr3d_trans mr3d_filter mr3d_recons mr3d_stat mr2d1d_trans) diff --git a/src/im1d_convert.cc b/src/im1d_convert.cc new file mode 100644 index 0000000..8f4e380 --- /dev/null +++ b/src/im1d_convert.cc @@ -0,0 +1,123 @@ +/****************************************************************************** +** Copyright (C) 1998 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 21/09/98 +** +** File: im1d_convert.cc +** +******************************************************************************* +** +** DESCRIPTION +** ----------- +** +** +** +******************************************************************************/ + + +#include "Array.h" +#include "IM_Obj.h" +#include "IM1D_IO.h" + +char Name_Imag_In[256]; +char Name_Imag_Out[256]; + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char *const*argv, char *opts); + +int NSkip=0; + +/****************************************************************************/ + +static void usage(char *argv[]) +{ + + fprintf(stdout, "Usage: %s options image output\n\n", argv[0]); + fprintf(stdout, " FITS format (.fits) \n"); + fprintf(stdout, " ACII format (.dat) \n"); + fprintf(stdout, " EXEL format (.csv) \n"); + + fprintf(stdout, " where options = \n"); + + fprintf(stdout, " [-s line_number_to_skip]\n"); + fprintf(stdout, " Skip the s first lines of the input file.\n"); + + fprintf(stdout, "\n"); + + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ + +static void transinit(int argc, char *argv[]) +{ + int c; + + /* get options */ + while ((c = GetOpt(argc,argv,"s:")) != -1) + { + switch (c) + { + case 's': + if (sscanf(OptArg,"%d",&NSkip) != 1) + { + fprintf(stdout, "bad number of lines: %s\n", OptArg); + usage(argv); + } + break; + case '?': + usage(argv); + } + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(stdout, "Too many parameters: %s ...\n", argv[OptInd]); + usage(argv); + } +} + + +/****************************************************************************/ + + +int main(int argc, char *argv[]) +{ + fltarray Dat; + fltarray Mallat; + /* Get command line arguments, open input file(s) if necessary */ + transinit(argc, argv); + + // Dat.fits_read(Name_Imag_In); + io_1d_read_data (Name_Imag_In, Dat, NSkip); + cout << endl << endl; + + cerr << "CONVERT " << String1DFormat(io_detect_1dformat(Name_Imag_In)); + cerr << " to " << String1DFormat(io_detect_1dformat(Name_Imag_Out)) << endl; + + io_1d_set_format(Name_Imag_Out); + io_1d_write_data (Name_Imag_Out, Dat); + + exit(0); +} + diff --git a/src/im1d_der_loggamma.cc b/src/im1d_der_loggamma.cc new file mode 100644 index 0000000..13884f2 --- /dev/null +++ b/src/im1d_der_loggamma.cc @@ -0,0 +1,608 @@ +/****************************************************************************** +** Copyright (C) 1999 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 10/02/00 +** +** File: dct.cc +** +******************************************************************************* +** +** DESCRIPTION dct program +** ----------- +** +******************************************************************************/ + +#include "GlobalInc.h" +#include "IM_IO.h" + +char Name_Imag_In[256]; /* input file image */ +char Name_Imag_Out[256]; /* output file name */ + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char **argv, char *opts); + +Bool Verbose=False; +int Order=0; + +/***************************************/ + +static void usage(char *argv[]) +{ + // int i; + fprintf(OUTMAN, "Usage: %s options in_catalogue result\n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + + fprintf(OUTMAN, " [-o Order]\n"); + fprintf(OUTMAN, " Order value.\n"); + fprintf(OUTMAN, " Default is 0.\n"); + manline(); + + fprintf(OUTMAN, " [-v]\n"); + fprintf(OUTMAN, " Verbose.\n"); + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ +static void filtinit(int argc, char *argv[]) +{ + int c; + + /* get options */ + while ((c = GetOpt(argc,argv,"o:vzZ")) != -1) + { + switch (c) + { + case 'o': + if (sscanf(OptArg,"%d",&Order) != 1) + { + fprintf(OUTMAN, "Error: bad order value: %s\n", OptArg); + exit(-1); + } + break; + case 'v': Verbose = True;break; + case '?': usage(argv); break; + default: usage(argv); break; + } + } + + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } +} + + +/*******************************************************************/ +/*******************************************************************/ + + +static double s_PolyGamma(double x, int order); + +/*********************************************************************/ +/** @file polygamma.c + * Definitions for portable math library + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#define LOGDERIV_ORDER_MAX 4 +#define POLYGAMMA_ORDER_MAX LOGDERIV_ORDER_MAX +#define DBL_EPSILON 2.2204460492503131e-16 +#define NCBIMATH_PI 3.1415926535897932384626433832795 +#define NCBIMATH_LNPI 1.1447298858494001741434273513531 +#define NCBIMATH_LN2 0.69314718055994530941723212145818 +#define DIM(A) (sizeof(A)/sizeof((A)[0])) +#define ABS(a) ((a)>=0?(a):-(a)) +// #define MAX(a,b) ((a)>=(b)?(a):(b)) +// #define MIN(a,b) ((a)<=(b)?(a):(b)) + +/** Compute ln(abs(gamma(x))) to 10-digit accuracy + * @param x Point to evaluate ln(abs(gamma(x))) + * @return The function value + */ +static double s_LnGamma(double x) +{ + return lgamma(x); +} +/** Tabulated values of the first few factorials */ +static const double kPrecomputedFactorial[] = { + 1., 1., 2., 6., 24., 120., 720., 5040., 40320., 362880., 3628800., + 39916800., 479001600., 6227020800., 87178291200., 1307674368000., + 20922789888000., 355687428096000., 6402373705728000., + 121645100408832000., 2432902008176640000., 51090942171709440000., + 1124000727777607680000., 25852016738884976640000., + 620448401733239439360000., 15511210043330985984000000., + 403291461126605635584000000., 10888869450418352160768000000., + 304888344611713860501504000000., 8841761993739701954543616000000., + 265252859812191058636308480000000., 8222838654177922817725562880000000., + 263130836933693530167218012160000000., + 8683317618811886495518194401280000000., + 295232799039604140847618609643520000000. + }; + +double BLAST_Factorial(int n) +{ + if (n < 0) + return 0.0; /* Undefined! */ + + if (n < DIM(kPrecomputedFactorial)) + return kPrecomputedFactorial[n]; + + return exp(s_LnGamma((double)(n + 1))); +} + +double BLAST_LnGammaInt(int n) +{ + if ( (n > 1) && (n < DIM(kPrecomputedFactorial) ) ) { + return log(kPrecomputedFactorial[n-1]); + } + + return s_LnGamma((double)n); +} + +/* + Romberg numerical integrator + Reference: + Francis Scheid (1968) + Schaum's Outline Series + Numerical Analysis, p. 115 + McGraw-Hill Book Company, New York +*/ + +/** Make a parametrized function appear to have only one variable */ +#define F(x) ((*f)((x), fargs)) +/** Maximum number of diagonals in the Romberg array */ +#define MAX_DIAGS 20 + +double BLAST_RombergIntegrate(double (*f) (double,void*), void* fargs, double p, double q, double eps, int epsit, int itmin) + +{ + double romb[MAX_DIAGS]; /* present list of Romberg values */ + double h; /* mesh-size */ + int i, j, k, npts; + long n; /* 4^(error order in romb[i]) */ + int epsit_cnt = 0, epsck; + double y; + double x; + double sum; + + /* itmin = min. no. of iterations to perform */ + itmin = MAX(1, itmin); + itmin = MIN(itmin, MAX_DIAGS-1); + + /* epsit = min. no. of consecutive iterations that must satisfy epsilon */ + epsit = MAX(epsit, 1); /* default = 1 */ + epsit = MIN(epsit, 3); /* if > 3, the problem needs more prior analysis */ + + epsck = itmin - epsit; + + npts = 1; + h = q - p; + x = F(p); + if (ABS(x) == HUGE_VALF) + return x; + y = F(q); + if (ABS(y) == HUGE_VALF) + return y; + romb[0] = 0.5 * h * (x + y); /* trapezoidal rule */ + for (i = 1; i < MAX_DIAGS; ++i, npts *= 2, h *= 0.5) { + sum = 0.; /* sum of ordinates for + x = p+0.5*h, p+1.5*h, ..., q-0.5*h */ + for (k = 0, x = p+0.5*h; k < npts; k++, x += h) { + y = F(x); + if (ABS(y) == HUGE_VALF) + return y; + sum += y; + } + romb[i] = 0.5 * (romb[i-1] + h*sum); /* new trapezoidal estimate */ + + /* Update Romberg array with new column */ + for (n = 4, j = i-1; j >= 0; n *= 4, --j) + romb[j] = (n*romb[j+1] - romb[j]) / (n-1); + + if (i > epsck) { + if (ABS(romb[1] - romb[0]) > eps * ABS(romb[0])) { + epsit_cnt = 0; + continue; + } + ++epsit_cnt; + if (i >= itmin && epsit_cnt >= epsit) + return romb[0]; + } + } + return HUGE_VALF; +} + +int BLAST_Gcd(int a, int b) +{ + int c; + + b = ABS(b); + if (b > a) + c=a, a=b, b=c; + + while (b != 0) { + c = a%b; + a = b; + b = c; + } + return a; +} + +int +BLAST_Gdb3(int* a, int* b, int* c) +{ + int g; + if (*b == 0) + g = BLAST_Gcd(*a, *c); + else + g = BLAST_Gcd(*a, BLAST_Gcd(*b, *c)); + if (g > 1) { + *a /= g; + *b /= g; + *c /= g; + } + return g; +} + +long BLAST_Nint(double x) +{ + x += (x >= 0. ? 0.5 : -0.5); + return (long)x; +} + + +double BLAST_Powi(double x, int n) +{ + double y; + + if (n == 0) + return 1.; + + if (x == 0.) { + if (n < 0) { + return HUGE_VALF; + } + return 0.; + } + + if (n < 0) { + x = 1./x; + n = -n; + } + + y = 1.; + while (n > 0) { + if (n & 1) + y *= x; + n /= 2; + x *= x; + } + return y; + +} + +double BLAST_LnFactorial (double x) { + + if( x <= 0.0) + return 0.0; + else + return s_LnGamma(x + 1.0); + +} + +double BLAST_Expm1(double x) +{ + double absx = ABS(x); + + if (absx > .33) + return exp(x) - 1.; + + if (absx < 1.e-16) + return x; + + return x * (1. + x * + (1./2. + x * + (1./6. + x * + (1./24. + x * + (1./120. + x * + (1./720. + x * + (1./5040. + x * + (1./40320. + x * + (1./362880. + x * + (1./3628800. + x * + (1./39916800. + x * + (1./479001600. + + x/6227020800.)))))))))))); +} + +/** size of the next series term that indicates convergence + in the log and polygamma functions */ +#ifndef DBL_EPSILON +#define DBL_EPSILON 2.2204460492503131e-16 +#endif + +double BLAST_Log1p(double x) +{ + int i; + double sum, y; + + if (ABS(x) >= 0.2) + return log(x+1.); + + /* Limit the loop to 500 terms. */ + for (i=0, sum=0., y=x; i<500 ; ) { + sum += y/++i; + if (ABS(y) < DBL_EPSILON) + break; + y *= x; + sum -= y/++i; + if (y < DBL_EPSILON) + break; + y *= x; + } + return sum; +} + +/** evaluate a specified-order derivative of ln(f(x)) + * @param order The order of derivative to evaluate (0...LOGDERIV_ORDER_MAX). + * Derivative order 0 just computes ln(f(x)) + * @param u A list of numerical values of f(x) and its derivatives, all at + * the same point x, to be used within the computations + * @return 'order'-th derivative of ln(f(x)) or HUGE_VALF if + * order is out of range or u[0] is zero + */ +static double s_LogDerivative(int order, double* u) +{ + int i; + double y[LOGDERIV_ORDER_MAX+1]; + double value, tmp; + + if (order < 0 || order > LOGDERIV_ORDER_MAX) { + return HUGE_VALF; + } + + if (order > 0 && u[0] == 0.) { + return HUGE_VALF; + } + for (i = 1; i <= order; i++) + y[i] = u[i] / u[0]; + + switch (order) { + case 0: + if (u[0] > 0.) + value = log(u[0]); + else { + return HUGE_VALF; + } + break; + case 1: + value = y[1]; + break; + case 2: + value = y[2] - y[1] * y[1]; + break; + case 3: + value = y[3] - 3. * y[2] * y[1] + 2. * y[1] * y[1] * y[1]; + break; + case 4: + value = y[4] - 4. * y[3] * y[1] - 3. * y[2] * y[2] + + 12. * y[2] * (tmp = y[1] * y[1]); + value -= 6. * tmp * tmp; + break; + default: + return HUGE_VALF; + } + return value; +} + +/** auxiliary values for computation of derivative of ln(gamma(x)) */ +static double _default_gamma_coef [] = { + 4.694580336184385e+04, + -1.560605207784446e+05, + 2.065049568014106e+05, + -1.388934775095388e+05, + 5.031796415085709e+04, + -9.601592329182778e+03, + 8.785855930895250e+02, + -3.155153906098611e+01, + 2.908143421162229e-01, + -2.319827630494973e-04, + 1.251639670050933e-10 + }; + +/** Compute a specified-order derivative of ln(gamma(x)) + * evaluated at some point x + * @param x Value at which derivative will be evaluated + * @param order Order of derivative (0...POLYGAMMA_ORDER_MAX) + * @return 'order'-th derivative of ln(gamma(x)) at specified x. + * Accuracy is to 10 digits for x >= 1 + */ +static double +s_GeneralLnGamma(double x, int order) +{ + int i; + double xx, tx; + double y[POLYGAMMA_ORDER_MAX+1]; + double tmp, value; + double *coef; + const int xgamma_dim = DIM(_default_gamma_coef); + + + xx = x - 1.; /* normalize from gamma(x + 1) to xx! */ + tx = xx + xgamma_dim; + for (i = 0; i <= order; ++i) { + tmp = tx; + /* sum the least significant terms first */ + coef = &_default_gamma_coef[xgamma_dim]; + if (i == 0) { + value = *--coef / tmp; + while (coef > _default_gamma_coef) + value += *--coef / --tmp; + } + else { + value = *--coef / BLAST_Powi(tmp, i + 1); + while (coef > _default_gamma_coef) + value += *--coef / BLAST_Powi(--tmp, i + 1); + tmp = BLAST_Factorial(i); + value *= (i%2 == 0 ? tmp : -tmp); + } + y[i] = value; + } + ++y[0]; + value = s_LogDerivative(order, y); + tmp = tx + 0.5; + switch (order) { + case 0: + value += ((NCBIMATH_LNPI+NCBIMATH_LN2) / 2.) + + (xx + 0.5) * log(tmp) - tmp; + break; + case 1: + value += log(tmp) - xgamma_dim / tmp; + break; + case 2: + value += (tmp + xgamma_dim) / (tmp * tmp); + break; + case 3: + value -= (1. + 2.*xgamma_dim / tmp) / (tmp * tmp); + break; + case 4: + value += 2. * (1. + 3.*xgamma_dim / tmp) / (tmp * tmp * tmp); + break; + default: + tmp = BLAST_Factorial(order - 2) * BLAST_Powi(tmp, 1 - order) + * (1. + (order - 1) * xgamma_dim / tmp); + if (order % 2 == 0) + value += tmp; + else + value -= tmp; + break; + } + return value; +} + +/** Compute, to 10-digit accuracy, a specified order + * derivative of ln(abs(gamma(x))). + * @param x value at which derivative will be evaluated + * @param order Order of derivative (0...POLYGAMMA_ORDER_MAX) + * order = 0, 1, 2, corresponds to ln(gamma), + * digamma, trigamma, etc. Note that the value here + * is one less than that suggested by the "di" and "tri" + * prefixes of digamma, trigamma, etc. In other words, + * it is truly the order of the derivative. + * @return Computed derivative value, or HUGE_VALF if order is out of range + */ +static double s_PolyGamma(double x, int order) +{ + int k; + double value, tmp; + double y[POLYGAMMA_ORDER_MAX+1], sx; + + order+=1; // To match the notation of Abramovitz and Stegun. + + if (order < 0 || order > POLYGAMMA_ORDER_MAX) { + return HUGE_VALF; + } + + if (x >= 1.) + return s_GeneralLnGamma(x, order); + + if (x < 0.) { + value = s_GeneralLnGamma(1. - x, order); + value = ((order - 1) % 2 == 0 ? value : -value); + if (order == 0) { + sx = sin(NCBIMATH_PI * x); + sx = ABS(sx); + if ( (x < -0.1 && (ceil(x) == x || sx < 2.*DBL_EPSILON)) + || sx == 0.) { + return HUGE_VALF; + } + value += NCBIMATH_LNPI - log(sx); + } + else { + y[0] = sin(x *= NCBIMATH_PI); + tmp = 1.; + for (k = 1; k <= order; k++) { + tmp *= NCBIMATH_PI; + y[k] = tmp * sin(x += (NCBIMATH_PI/2.)); + } + value -= s_LogDerivative(order, y); + } + } + else { + value = s_GeneralLnGamma(1. + x, order); + if (order == 0) { + if (x == 0.) { + return HUGE_VALF; + } + value -= log(x); + } + else { + tmp = BLAST_Factorial(order - 1) * BLAST_Powi(x, -order); + value += (order % 2 == 0 ? tmp : - tmp); + } + } + return value; +} + + + + +int main(int argc, char *argv[]) +{ + int k; + fitsstruct Header; + char Cmd[512]; + Cmd[0] = '\0'; + for (k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + + /* Get command line arguments, open input file(s) if necessary */ + filtinit(argc, argv); + + if (Verbose == True) + { + cout << endl << endl << "PARAMETERS: " << endl << endl; + cout << "File Name in = " << Name_Imag_In << endl; + cout << "File Name Out = " << Name_Imag_Out << endl; + cout << "Order = " << Order << endl; + } + + fltarray Data; + fits_read_fltarr(Name_Imag_In, Data, &Header); + for (int i=0; i < Data.nx(); i++) Data(i) = s_PolyGamma(Data(i), Order); + fits_write_fltarr(Name_Imag_Out, Data); + exit(0); +} diff --git a/src/im1d_info.cc b/src/im1d_info.cc new file mode 100644 index 0000000..1b9233e --- /dev/null +++ b/src/im1d_info.cc @@ -0,0 +1,201 @@ +/****************************************************************************** +** Copyright (C) 2000 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 28/08/00 +** +** File: im1d_info.cc +** +******************************************************************************* +** +** DESCRIPTION return information about a time serie data set. +** ----------- +** +******************************************************************************/ + +#include "Array.h" +#include "IM_Obj.h" +#include "IM_IO.h" +#include "IM1D_IO.h" +#include "Usage.h" +#include "MatrixOper.h" +#ifdef USELM +#include "Licence.h" +#endif + +char Name_Imag_In[256]; +char Name_Imag_Out[256]; +char Name_Imag_Out2[256]; +char Name_Coeff[256]; + + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char **argv, char *opts); + +Bool Verbose = False; +int NBShift = 10; +ar_order_detect_type AR_Order_Detect=AR_ORDER_BIC; +int BestNbrAR; +int MaxNbrAR=10; +Bool AutoCorr = False; + +/***************************************************************************/ + +static void usage(char *argv[]) +{ + + fprintf(OUTMAN, "Usage: %s options input \n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + + fprintf(OUTMAN, "\n"); + fprintf(OUTMAN, " [-a Nbr_of_Lags]\n"); + fprintf(OUTMAN, " Calculate the autocorrelation function (AF)\n"); + fprintf(OUTMAN, " with a given number of lags.\n"); + fprintf(OUTMAN, " The output AF is saved in a file.\n"); + fprintf(OUTMAN, " File Name is: autocor.\n"); + // fprintf(OUTMAN, " Default number of lags is %d.\n", NBShift); + manline(); + + fprintf(OUTMAN, " [-O Estimation_AR_Order_Method]\n"); + fprintf(OUTMAN, " 1: AIC \n"); + fprintf(OUTMAN, " 2: AICC\n"); + fprintf(OUTMAN, " 3: BIC\n"); + fprintf(OUTMAN, " Default is %s method.\n", StringARDetectType(AR_Order_Detect)); + manline(); + + fprintf(OUTMAN, " [-M MaxAROrder]\n"); + fprintf(OUTMAN, " Default is %d.\n", MaxNbrAR); + manline(); + verbose_usage(); + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ +static void filtinit(int argc, char *argv[]) +{ + int c; + /* get options */ + while ((c = GetOpt(argc,argv,"a:O:M:v")) != -1) + { + switch (c) + { + case 'O': + if (sscanf(OptArg,"%d",&c) != 1) + { + fprintf(OUTMAN, "Error: bad order method: %s\n", OptArg); + exit(-1); + } + AR_Order_Detect = (ar_order_detect_type)(c-1); + break; + case 'M': + if (sscanf(OptArg,"%d",&MaxNbrAR) != 1) + { + fprintf(OUTMAN, "Error: bad maximum order value: %s\n", OptArg); + exit(-1); + } + break; + case 'a': + if (sscanf(OptArg,"%d",&NBShift) != 1) + { + fprintf(OUTMAN, "Error: bad number: %s\n", OptArg); + exit(-1); + } + AutoCorr = True; + strcpy(Name_Imag_Out, "autocor"); + break; + case 'v': Verbose = True; break; + case '?': usage(argv); break; + } + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Too many parameters: %s ...\n", argv[OptInd]); + usage(argv); + } +} + +/****************************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Data,TrueData; + char Cmd[256]; + // extern softinfo Soft; + +#ifdef USELM + lm_check(LIC_M1D); +#endif + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + /* Get command line arguments, open input file(s) if necessary */ + + filtinit(argc, argv); + + io_1d_read_data(Name_Imag_In, Data); + // reform_to_1d(Data); + int Nx = Data.nx(); + if (NBShift >= Nx-10) + { + cout << "Error: NBShift = " << NBShift << endl; + cout << " 1 < NBShift < " << Nx -10 << endl; + } + + // FitsHeader.origin = Cmd; + + if (Verbose == True) + { + cout << endl << endl << "PARAMETERS: " << endl << endl; + cout << "File Name in = " << Name_Imag_In << endl; + if (AutoCorr == True) + cout << "File Name out = " << Name_Imag_Out << endl; + cout << "AR Order Detection Method = " << StringARDetectType(AR_Order_Detect) << endl; + if (AutoCorr == True) + cout << "AF calculation: NBLags = " << NBShift << endl << endl; + } + + cout << " Np = " << Data.nx() << endl; + float FluxIma = Data.total(); + double Mean,Sigma,Skew,Curt; + float Min,Max; + moment4(Data.buffer(),Data.nx(), Mean, Sigma, Skew, Curt, Min,Max); + cout << " Min = " << Min << " Max = " << Max << endl; + cout << " Mean = " << Mean << " sigma = " << Sigma << endl; + cout << " Skew = " << Skew << " Curtosis = " << Curt << endl; + cout << " Flux = " << FluxIma << " Energy = " << (Data^2).total() << endl; + + fltarray BestARModel; + AR_PREDICTION AR; + AR.AR_Order_Detect = AR_Order_Detect; + AR.MaxNbrAR = MaxNbrAR; + AR.Verbose = Verbose; + + if (AutoCorr == True) + { + fltarray TabAuto(NBShift); + autocor1d(Data, TabAuto, NBShift); + io_1d_write_data(Name_Imag_Out, TabAuto); + } + AR.get_best_ar_model(Data, Data.nx(), BestNbrAR, BestARModel); + cout << " Optimal AR order = " << BestNbrAR << endl; + + exit(0); +} + diff --git a/src/im1d_period.cc b/src/im1d_period.cc new file mode 100644 index 0000000..12f25d3 --- /dev/null +++ b/src/im1d_period.cc @@ -0,0 +1,398 @@ +/****************************************************************************** +** Copyright (C) 2003 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Philippe Querre & Jean-Luc Starck +** +** Date: 19/03/03 +** +** File: im1d_period.cc +** +******************************************************************************* +** +** DESCRIPTION periodogram calcualation +** ----------- +** +******************************************************************************/ + + + +#include "IM_Obj.h" +#include "IM_IO.h" +#include "IM1D_IO.h" +#include "FFTN_1D.h" + +extern int GetOpt(int argc, char **argv, char *opts); +extern int OptInd; +extern char *OptArg; + + +char Name_Imag_In[256]; /* input file image */ +char Name_Imag_Out[256]; /* output file name */ +Bool Verbose=False; +Bool SubstractMean = False; +int DanielHalfSize=5; +int BarlettSegmentNumber=8; +int WelchShift=0; +int WelchNbEchant=0; +int NbPointOut=0; +float WinParam = 0.5; +type_std_win WinType = W_HANNING; + +#define NBR_PERIODOGRAM_TYPE 4 +enum periodogram_type {STANDARD=1, + BARLETT=2, + WELCH=3, + DANIEL=4}; +periodogram_type PeriodogramType=WELCH; + +inline char * StringPeriodogramType (periodogram_type PT) +{ + switch (PT) + { + case STANDARD : return ("Standard periodogram");break; + case BARLETT : return ("Barlett periodogram");break; + case WELCH : return ("Welch periodogram");break; + case DANIEL : return ("Daniel periodogram");break; + default: return ("Undefined periodogram type"); + exit(-1);break; + } +} + + +/***************************************/ + +static void usage(char *argv[]) { + + + + fprintf(OUTMAN, "Usage: %s options Reference Data\n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + manline(); + + fprintf(OUTMAN, " [-m Periodogram method]\n"); + for (int i=1;i<=NBR_PERIODOGRAM_TYPE;i++) + fprintf(OUTMAN, " %d: %s \n",i,StringPeriodogramType((periodogram_type)i)); + fprintf(OUTMAN, " default is 3.\n"); + manline(); + + fprintf(OUTMAN, " [-h HalfSizeWindow]\n"); + fprintf(OUTMAN, " Window half size used in the Daniel periodogram.\n"); + fprintf(OUTMAN, " default is 5.\n"); + manline(); + + fprintf(OUTMAN, " [-q SegmentNbr]\n"); + fprintf(OUTMAN, " Number of segment used in the Barlett periodogram.\n"); + fprintf(OUTMAN, " default is 8\n"); + manline(); + + fprintf(OUTMAN, " [-r Shift]\n"); + fprintf(OUTMAN, " Shift used in Welch periodogram.\n"); + fprintf(OUTMAN, " default is signal_point_number / 10\n"); + manline(); + + fprintf(OUTMAN, " [-s SegmentNbr]\n"); + fprintf(OUTMAN, " Number of segment used in Welch periodogram.\n"); + fprintf(OUTMAN, " default is signal_point_number / 2\n"); + manline(); + + fprintf(OUTMAN, " [-t PeriodPointNbr]\n"); + fprintf(OUTMAN, " Number of point of periodogram.\n"); + fprintf(OUTMAN, " It must be greater or equal to signal_point_number\n"); + fprintf(OUTMAN, " default is signal_point_number.\n"); + manline(); + + fprintf(OUTMAN, " [-T type_of_window ]\n"); + for (int i = 0; i < NBR_STD_WIN ; i++) + fprintf(OUTMAN, " %d: %s \n",i+1,StringSTDWin((type_std_win )i)); + fprintf(OUTMAN, " default is %s.\n", StringSTDWin(WinType)); + manline(); + + fprintf(OUTMAN, " [-W window_param (time domain)]\n"); + fprintf(OUTMAN, " Window parameter. Default is 0.5.\n"); + manline(); + + fprintf(OUTMAN, " [-M]"); + fprintf(OUTMAN, " Substract mean (default is no)"); + manline(); + + verbose_usage(); + + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ +static void infinit(int argc, char *argv[]) { + + int c; +#ifdef LARGE_BUFF + int VMSSize=-1; + Bool OptZ = False; + char VMSName[1024] = ""; +#endif + + /* get options */ + int v; + while ((c = GetOpt(argc,argv,"m:h:q:r:s:t:T:W:MvzZ")) != -1) { + + switch (c) { + case 'v': Verbose = True; break; + case 'M': SubstractMean = True; break; + case 'm': + if (sscanf(OptArg,"%d",&v ) != 1) { + fprintf(OUTMAN, "bad periodogram type !!! : %s\n", OptArg); + exit(-1); + } + if ((v<1) || (v>NBR_PERIODOGRAM_TYPE)) + PeriodogramType = WELCH; + else PeriodogramType = (periodogram_type)v; + cout << StringPeriodogramType(PeriodogramType) << endl; + break; + case 'h': + if (sscanf(OptArg,"%d",&DanielHalfSize) != 1) { + fprintf(OUTMAN, "bad Daniel periodogram half size value!!! : %s\n", + OptArg); + exit(-1); + } + break; + case 'q': + if (sscanf(OptArg,"%d",&BarlettSegmentNumber) != 1) { + fprintf(OUTMAN, "bad Barlett segment number!!! : %s\n", + OptArg); + exit(-1); + } + break; + case 'r': + if (sscanf(OptArg,"%d",&WelchShift) != 1) { + fprintf(OUTMAN, "bad Welch shift!!! : %s\n", + OptArg); + exit(-1); + } + break; + case 's': + if (sscanf(OptArg,"%d",&WelchNbEchant) != 1) { + fprintf(OUTMAN, "bad Welch point number!!! : %s\n", + OptArg); + exit(-1); + } + break; + case 't': + if (sscanf(OptArg,"%d",&NbPointOut) != 1) { + fprintf(OUTMAN, "bad out point number!!! : %s\n", + OptArg); + exit(-1); + } + break; + case 'T': + if (sscanf(OptArg,"%d",&v) != 1) { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + exit(-1); + } + if ((v <= 0) && (v > NBR_STD_WIN)) { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + usage(argv); + } + WinType = (type_std_win) (v-1); + break; + case 'W': + if (sscanf(OptArg,"%f",&WinParam) != 1) { + fprintf(OUTMAN, "bad window parameter: %s\n", OptArg); + exit(-1); + } + break; + +#ifdef LARGE_BUFF + case 'z': + if (OptZ == True) + { + fprintf(OUTMAN, "Error: Z option already set...\n"); + exit(-1); + } + OptZ = True; + break; + case 'Z': + if (sscanf(OptArg,"%d:%s",&VMSSize, VMSName) < 1) + { + fprintf(OUTMAN, "Error: syntaxe is Size:Directory ... \n"); + exit(-1); + } + if (OptZ == True) + { + fprintf(OUTMAN, "Error: z option already set...\n"); + exit(-1); + } + OptZ = True; + break; +#endif + case '?': usage(argv); break; + default: usage(argv); break; + } + } + + + /* get optional input file names from trailing parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) { + fprintf(OUTMAN, "Too many parameters: %s ...\n", argv[OptInd]); + exit(-1); + } + + +#ifdef LARGE_BUFF + if (OptZ == True) vms_init(VMSSize, VMSName, Verbose); +#endif +} + + +/****************************************************************************/ + +int main(int argc, char *argv[]) +{ + /* Get command line arguments, open input file(s) if necessary */ + lm_check(LIC_MR1); + infinit(argc, argv); + Ifloat IDat; + fltarray Data; + + io_1d_read_data(Name_Imag_In, Data); + int Nx = Data.nx(); + if (SubstractMean) { + float Mean = Data.mean(); + for (int i=0;i NbEchant) // add sero (zero padding) + for (int j=NbEchant;j NBR_STD_WIN)) + { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + usage(argv); + } + WindowType = (type_std_win) (c-1); + break; + case 'W': + /* -n */ + if (sscanf(OptArg,"%f",&WinParam) != 1) + { + fprintf(OUTMAN, "bad window parameter: %s\n", OptArg); + exit(-1); + } + break; + case 'w': + /* -n */ + if (sscanf(OptArg,"%d",&WindowSize) != 1) + { + fprintf(OUTMAN, "bad window size: %s\n", OptArg); + exit(-1); + } + break; + case 'S': + /* -n */ + if (sscanf(OptArg,"%d",&Step) != 1) + { + fprintf(OUTMAN, "bad step: %s\n", OptArg); + exit(-1); + } + break; + case 'r': Reverse = True; break; + case 'v': Verbose = True; break; + case '?': usage(argv); break; + } + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Too many parameters: %s ...\n", argv[OptInd]); + usage(argv); + } +} + + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Data; + Ifloat STF_Spec; + Icomplex_f STFDat; + char Cmd[256]; + int i,j,Nx,Nlt,Nct; + char Name_Imag_Out_Spec[256]; + ST_FFTN STF; + complex_f *Buff_STF; + + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + /* Get command line arguments, open input file(s) if necessary */ + + lm_check(LIC_M1D); + filtinit(argc, argv); + + if (Reverse == False) + { + io_1d_read_data(Name_Imag_In, Data); + // reform_to_1d(Data); + // FitsHeader.origin = Cmd; + Nx = Data.nx(); + if (WindowSize <= 0) WindowSize = Nx / 4; + else if (WindowSize >= Nx) WindowSize = Nx / 4; + if (Step <= 0) Step = WindowSize / 2; + STF.alloc(Nx, WindowType, WinParam, WindowSize, Step); + Nlt = STF.nl(Nx); + Nct = STF.nc(Nx); + STFDat.alloc(Nlt, Nct, "STF"); + } + else + { + //cout << "READ COMPLEX: " << Name_Imag_In << endl; + io_read_ima_complex_f (Name_Imag_In, STFDat); + //cout << "CF: " << STFDat.nl() << " " << STFDat.nc() << endl; + Nlt = STFDat.nl(); + Nct = STFDat.nc(); + WindowSize = Nlt; + if (Step <= 0) Step = WindowSize / 2; + Nx = Nct*Step; + //cout << "ALLOC" << Nx << " WindowSize = " << WindowSize << " Step = " << Step << endl; + STF.alloc(Nx, WindowType, WinParam, WindowSize, Step); + Data.alloc(Nx); + } + Buff_STF = STFDat.buffer(); + + if (Verbose == True) + { + cout << endl << endl << "PARAMETERS: " << endl << endl; + cout << "File Name in = " << Name_Imag_In << endl; + cout << "File Name out = " << Name_Imag_Out << endl; + cout << "Window Type = " << StringSTDWin(WindowType) << endl; + cout << "Window size = " << WindowSize << endl; + cout << "Step = " << Step << endl; + cout << "Window parameter = " << WinParam << endl; + cout << "Data size = " << Nx << endl; + cout << "STF size: Nl = " << Nlt << " Nc = " << Nct << endl; + } + + if (DebugSTF == True) + { + float W=WinParam; + fltarray Win(Nx); + STD_Window STDW; + STDW.hamming(Win, W); + io_1d_write_data("xx_hamming.fits", Win); + STDW.hanning(Win, W); + io_1d_write_data("xx_hanning.fits", Win); + STDW.gaussian(Win, W); + io_1d_write_data("xx_gaussian.fits", Win); + STDW.blackman(Win, W); + io_1d_write_data("xx_blackman.fits", Win); + } + + if (Reverse == False) + { + STF.transform(Data, Buff_STF); + io_write_ima_complex_f (Name_Imag_Out, STFDat); + STF_Spec.alloc(Nlt/2,Nct,"STF spec"); + for (i=0; i < Nlt/2; i++) + for (j=0; j < Nct; j++) + { + float Re = Buff_STF[i*Nct+j].real(); // / (float) Nx; + float Im = Buff_STF[i*Nct+j].imag(); // / (float) Nx; + STF_Spec(i,j) = Re*Re+Im*Im; + } + io_write_ima_complex_f (Name_Imag_Out, STFDat); + sprintf(Name_Imag_Out_Spec,"%s_spec.fits", Name_Imag_Out); + io_write_ima_float(Name_Imag_Out_Spec, STF_Spec); + } + else + { + STF.recons(Buff_STF, Data, True); + io_1d_write_data(Name_Imag_Out, Data); + } + exit(0); +} + diff --git a/src/im1d_tend.cc b/src/im1d_tend.cc new file mode 100644 index 0000000..de45862 --- /dev/null +++ b/src/im1d_tend.cc @@ -0,0 +1,152 @@ +/****************************************************************************** +** Copyright (C) 2000 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 28/08/00 +** +** File: tend_est.cc +** +******************************************************************************* +** +** DESCRIPTION tendency estimation +** ----------- +** +** +******************************************************************************/ + +#include "Array.h" +#include "IM_Obj.h" +#include "IM_IO.h" +#include "IM1D_IO.h" +#include "Usage.h" +#include "MatrixOper.h" + +char Name_Imag_In[256]; +char Name_Imag_Out[256]; +char Name_Imag_Out2[256]; +char Name_Coeff[256]; + + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char **argv, char *opts); + +Bool Verbose = False; +int NbrPixTendancy = 100; +int Nt=-1; + +/***************************************************************************/ + +static void usage(char *argv[]) +{ + fprintf(OUTMAN, "Usage: %s options input out_tend out_signal_no_tend \n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + + fprintf(OUTMAN, "\n"); + + fprintf(OUTMAN, " [-T WindowSize_for_tendency_estimation]\n"); + fprintf(OUTMAN, " Default is %d.\n", NbrPixTendancy); + manline(); + fprintf(OUTMAN, " [-f FirstPixels]\n"); + fprintf(OUTMAN, " Default is the input signal size.\n"); + manline(); + verbose_usage(); + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ +static void filtinit(int argc, char *argv[]) +{ + int c; + /* get options */ + while ((c = GetOpt(argc,argv,"T:f:v")) != -1) + { + switch (c) + { + case 'T': + if (sscanf(OptArg,"%d",&NbrPixTendancy) != 1) + { + fprintf(OUTMAN, "Error: bad number: %s\n", OptArg); + exit(-1); + } + break; + case 'f': + if (sscanf(OptArg,"%d",&Nt) != 1) + { + fprintf(OUTMAN, "Error: bad number: %s\n", OptArg); + exit(-1); + } + break; + case 'v': Verbose = True; break; + case '?': usage(argv); break; + } + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out2, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Too many parameters: %s ...\n", argv[OptInd]); + usage(argv); + } +} + +/****************************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Data,TrueData; + char Cmd[256]; + // extern softinfo Soft; + + lm_check(LIC_M1D); + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + /* Get command line arguments, open input file(s) if necessary */ + + filtinit(argc, argv); + + io_1d_read_data(Name_Imag_In, Data); + // reform_to_1d(Data); + int Nx = Data.nx(); + + if (Nt < 1) Nt = Nx; + else if (Nt > Nx) Nt = Nx; + // FitsHeader.origin = Cmd; + + if (Verbose == True) + { + cout << endl << endl << "PARAMETERS: " << endl << endl; + cout << "File Name in = " << Name_Imag_In << endl; + cout << " NbrPixTendancy = " << NbrPixTendancy << endl; + } + + fltarray Tend(Nt); + tendancy_est(Data, Tend, NbrPixTendancy, Nt); + fltarray ST(Nt); + for (int i = 0; i < Nt; i++) ST(i) = Data(i) - Tend(i); + io_1d_write_data(Name_Imag_Out, Tend); + io_1d_write_data(Name_Imag_Out2, ST); + exit(0); +} + diff --git a/src/im1d_tfreq.cc b/src/im1d_tfreq.cc new file mode 100644 index 0000000..d20ae3a --- /dev/null +++ b/src/im1d_tfreq.cc @@ -0,0 +1,252 @@ +/****************************************************************************** +** Copyright (C) 2001 by CEA +******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Jean-Luc Starck +** +** Date: 13/04/01 +** +** File: im1d_stf.cc +** +******************************************************************************* +** +** DESCRIPTION Short term Fourier transform +** ----------- +** +** +******************************************************************************/ + + +#include "IM_Obj.h" +#include "IM_IO.h" +#include "IM_Noise.h" +#include "IM1D_IO.h" +#include "FFTN_1D.h" + +char Name_Imag_In[256]; +char Name_Imag_Out[256]; +char Name_Coeff[256]; + +extern int OptInd; +extern char *OptArg; + +extern int GetOpt(int argc, char **argv, char *opts); + +Bool Verbose = False; +Bool Reverse = False; +float WinParam = 0.5; +int Step = 0; +int WindowSize = 1024; +type_std_win WindowType = W_HAMMING; +Bool DebugSTF = False; + +type_spec TSpec = DEF_SPEC; +type_std_win WindowTypeFreq = W_HAMMING; + +/***************************************************************************/ + + +static void usage(char *argv[]) +{ + // int i; + + fprintf(OUTMAN, "Usage: %s options input \n\n", argv[0]); + fprintf(OUTMAN, " where options = \n"); + + spec_usage(); + manline(); + + fprintf(OUTMAN, " [-T type_of_window (time domain)]\n"); + for (int i = 0; i < NBR_STD_WIN ; i++) + fprintf(OUTMAN, " %d: %s \n",i+1, + StringSTDWin((type_std_win )i)); + fprintf(OUTMAN, " default is %s.\n", StringSTDWin(WindowType)); + manline(); + + // fprintf(OUTMAN, " [-F type_of_window (freq. domain)]\n"); +// fprintf(OUTMAN, " Only for Choi-Williams distribution.\n"); +// fprintf(OUTMAN, " default is %s.\n", StringSTDWin(WindowTypeFreq)); +// manline(); + + fprintf(OUTMAN, " [-w window_size (time domain)]\n"); + fprintf(OUTMAN, " Window size. Default is 1024.\n"); + manline(); + + fprintf(OUTMAN, " [-W window_param (time domain)]\n"); + fprintf(OUTMAN, " Window parameter. Default is 0.5.\n"); + manline(); + + fprintf(OUTMAN, " [-S Step]\n"); + fprintf(OUTMAN, " Step between two points. Default is WindowSize/2.\n"); + manline(); + + verbose_usage(); + manline(); + exit(-1); +} + +/*********************************************************************/ + +/* GET COMMAND LINE ARGUMENTS */ +static void filtinit(int argc, char *argv[]) +{ + int c; + /* get options */ + while ((c = GetOpt(argc,argv,"t:F:T:W:w:S:v")) != -1) + { + switch (c) + { + case 't': + if (sscanf(OptArg,"%d",&c) != 1) + { + fprintf(OUTMAN, "bad distribution type: %s\n", OptArg); + exit(-1); + } + if ((c <= 0) && (c > NBR_SPEC)) + { + fprintf(OUTMAN, "bad distribution type: %s\n", OptArg); + usage(argv); + } + TSpec = (type_spec ) (c-1); + break; + case 'T': + if (sscanf(OptArg,"%d",&c) != 1) + { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + exit(-1); + } + if ((c <= 0) && (c > NBR_STD_WIN)) + { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + usage(argv); + } + WindowType = (type_std_win) (c-1); + break; + case 'F': + if (sscanf(OptArg,"%d",&c) != 1) + { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + exit(-1); + } + if ((c <= 0) && (c > NBR_STD_WIN)) + { + fprintf(OUTMAN, "bad window type: %s\n", OptArg); + usage(argv); + } + WindowTypeFreq = (type_std_win) (c-1); + break; + case 'W': + /* -n */ + if (sscanf(OptArg,"%f",&WinParam) != 1) + { + fprintf(OUTMAN, "bad window parameter: %s\n", OptArg); + exit(-1); + } + break; + case 'w': + /* -n */ + if (sscanf(OptArg,"%d",&WindowSize) != 1) + { + fprintf(OUTMAN, "bad window size: %s\n", OptArg); + exit(-1); + } + break; + case 'S': + /* -n */ + if (sscanf(OptArg,"%d",&Step) != 1) + { + fprintf(OUTMAN, "bad step: %s\n", OptArg); + exit(-1); + } + break; + case 'v': Verbose = True; break; + case '?': usage(argv); break; + } + } + + /* get optional input file names from trailing + parameters and open files */ + if (OptInd < argc) strcpy(Name_Imag_In, argv[OptInd++]); + else usage(argv); + + if (OptInd < argc) strcpy(Name_Imag_Out, argv[OptInd++]); + else usage(argv); + + /* make sure there are not too many parameters */ + if (OptInd < argc) + { + fprintf(OUTMAN, "Too many parameters: %s ...\n", argv[OptInd]); + usage(argv); + } +} + + +/*********************************************************************/ + +int main(int argc, char *argv[]) +{ + fltarray Data; + fltarray TimeFreq; + Icomplex_f STFDat; + char Cmd[256]; + int Nx,Nlt,Nct; + ST_FFTN STF; + + Cmd[0] = '\0'; + for (int k =0; k < argc; k++) sprintf(Cmd, "%s %s", Cmd, argv[k]); + /* Get command line arguments, open input file(s) if necessary */ + + lm_check(LIC_M1D); + filtinit(argc, argv); + + io_1d_read_data(Name_Imag_In, Data); + // reform_to_1d(Data); + // FitsHeader.origin = Cmd; + Nx = Data.nx(); + if (WindowSize <= 0) WindowSize = Nx / 4; + else if (WindowSize >= Nx) WindowSize = Nx / 4; + if (Step <= 0) Step = WindowSize / 2; + + STF.alloc(Nx, WindowType, WinParam, WindowSize, Step); + Nlt = STF.nl(Nx); + Nct = STF.nc(Nx); + + if (Verbose == True) + { + cout << endl << endl << "PARAMETERS: " << endl << endl; + cout << "File Name in = " << Name_Imag_In << endl; + cout << "File Name out = " << Name_Imag_Out << endl; + cout << "Distribution type = " << StringSPEC(TSpec) << endl; + if (TSpec != SPEC_WIGNER_VILLE) + { + cout << "Window Type = " << StringSTDWin(WindowType) << endl; + cout << "Window size = " << WindowSize << endl; + cout << "Step = " << Step << endl; + cout << "Window parameter = " << WinParam << endl; + } + cout << "Data size = " << Nx << endl; + cout << "Output distrib. size: Nl = " << Nlt/2 << " Nc = " << Nct << endl; + } + + switch(TSpec) + { + case SPEC_STF: + STF.spectogram(Data, TimeFreq); + break; + case SPEC_WIGNER_VILLE: + STF.wigner_wille(Data, TimeFreq); + break; + case SPEC_CHOI_WILLIAMS: + STF.choi_williams(Data, TimeFreq, WindowTypeFreq); + break; + } + // cout << "Write " << Name_Imag_Out << TimeFreq.nl() << " " << TimeFreq.nc() << endl; + io_1d_write_data(Name_Imag_Out, TimeFreq); + + exit(0); +} + diff --git a/src/libsparse1d/MR1D_Deconv.cc b/src/libsparse1d/MR1D_Deconv.cc new file mode 100755 index 0000000..9871ba3 --- /dev/null +++ b/src/libsparse1d/MR1D_Deconv.cc @@ -0,0 +1,1336 @@ +/******************************************************************************* +** +** UNIT +** +** Version: 1.0 +** +** Author: Philippe Querre +** +** Date: 03/03/03 +** +** File: MR1D_Deconv.cc +** +******************************************************************************* +** +** DESCRIPTION signal deconvolution using multiresolution +** ----------- +** +** +***************************************************************************/ + +#include "MR1D_Deconv.h" +#include "DefFunc.h" + + + +/**************************************************************************/ +/********** ****************************************/ +/********** smooth_mediane ****************************************/ +/********** ****************************************/ +/**************************************************************************/ +void smooth_mediane (const fltarray& Sig1, fltarray& Sig2, + type_border Type, int Step_trou, int Window_Size) { + + int Nx = Sig1.nl(); + int k,i,ind_fen; + int Step; + float *fenetre; + int Window2 = (Window_Size - 1) / 2; + int Size = Window_Size; + + if (Step_trou > 0) Step = (int)(pow((double)2., (double) Step_trou) + 0.5); + else Step = 1; + + fenetre = new float [Size]; + fenetre[0]=0.; + for (i = 0; i < Nx; i++) { + + ind_fen = 0; + for (k = i - Window2*Step; k <= i+Window2*Step; k += Step) + fenetre[ind_fen++] = Sig1(k , Type); + + // Imag2(i,j) = hmedian(fenetre, ind_fen); + if (Size == 9) Sig2(i) = opt_med9(fenetre); + else Sig2(i) = get_median(fenetre, Size); + } + delete[] fenetre; +} + + +/**************************************************************************/ +/********** *********************************/ +/********** detect_noise_from_med *********************************/ +/********** *********************************/ +/**************************************************************************/ +/* search for an estimation of sigma in noise with median transform */ +/**************************************************************************/ +extern float BadPixalVal; +extern Bool BadPixel; +float detect_noise_from_med (const fltarray& Signal) { + + double Noise; + int i; + int N = Signal.nx(); + Bool Average_Non_Null = True; + int Nit = 3; + + fltarray Buff(Signal.nx(), "Buff noise estimation"); + smooth_mediane (Signal, Buff, DEFAULT_BORDER, 0, 3); + for (i=0; i < Signal.nx(); i++) Buff(i) = Signal(i) - Buff(i); + // INFO(Buff, "diff"); + + Noise=get_sigma_clip (Buff.buffer(), N, Nit, + Average_Non_Null, BadPixel, BadPixalVal); + Noise /= 0.972463; + // printf("Calc sigma(%d) = %f\n", N, (float) Noise); + // if (Noise < FLOAT_EPSILON) Noise = detect_noise_from_mad(Buff); + // INFO(Buff, "buff"); + return ((float) Noise); +} + +/**************************************************************************/ +/********** *******************************************/ +/********** dec_pos_max *******************************************/ +/********** *******************************************/ +/**************************************************************************/ +/* find max and ind_max of Tab in (only positive if SearchPositiv=True) */ +/**************************************************************************/ +void dec_pos_max (fltarray& Tab, int &Ind_i, float &Val_Max, + Bool SearchPositiv=False); +void dec_pos_max (fltarray& Tab, int &Ind_i, float &Val_Max, + Bool SearchPositiv) { + + int i; + float Val_Abs_Max; + int Bande = 1; + int Nx = Tab.nx(); + + Val_Max = 0.; + Val_Abs_Max = 0.; + + for (i = Bande; i < Nx-Bande; i++) { + + float Val_Abs = (SearchPositiv == False) ? ABS(Tab(i)) : Tab(i); + if (Val_Abs > Val_Abs_Max) { + + Val_Abs_Max = Val_Abs; + Val_Max = Tab(i); + Ind_i = i; + } + } +} + + +/**************************************************************************/ +/********** ****************************************/ +/********** dec_center_psf ****************************************/ +/********** ****************************************/ +/**************************************************************************/ +/* create a center psf of size Nx (signal in size), then normalize it */ +/**************************************************************************/ +void dec_center_psf (fltarray& D_Beam, fltarray& Psf) { + + int Ind_i,i,Dep_i; + float Val_Max; + int Nx_Beam = D_Beam.nx(); + int Nx = Psf.nx(); + double Flux = 0.; + + dec_pos_max (D_Beam, Ind_i, Val_Max, True); + + Psf.init(); + for (i = 0; i < Nx_Beam; i++) { + + Dep_i = i - Ind_i + Nx / 2; + if ((Dep_i >= 0) && (Dep_i < Nx)) { + + Psf(Dep_i) = D_Beam(i); + Flux += Psf(Dep_i); + } + } + for (i = 0; i < Nx; i++) Psf (i) = (float) (Psf(i) / Flux); +} + + + + +/**************************************************************************/ +/********** ********************************************/ +/********** psf_extend ********************************************/ +/********** ********************************************/ +/**************************************************************************/ +/* create a psf of signal size (but not center on max) */ +/**************************************************************************/ +void psf_extend (const fltarray& Signal, fltarray& Signal_Out) { + + int Nx0 = Signal.nx(); + int Nx1 = Signal_Out.nx(); + int i1,i0,Depi; + double Flux = 0.; + + Depi = (Nx1 - Nx0) / 2; + + for (i1 = 0; i1 < Nx1; i1++) { + + i0 = i1 - Depi; + if ((i0 < 0) || (i0 >= Nx0)) + Signal_Out (i1) = 0.; + else { + Signal_Out (i1) = Signal (i0); + Flux += Signal_Out (i1); + } + } + for (i1 = 0; i1 < Nx1; i1++) + Signal_Out (i1) = (float) (Signal_Out (i1) / Flux); +} + + +/**************************************************************************/ +/********** ***********************************************/ +/********** psf_get ***********************************************/ +/********** ***********************************************/ +/**************************************************************************/ +/* compute the fft of the psf (signal size) in Psf_cf */ +/**************************************************************************/ +void psf_get (fltarray& InPsf, cfarray &Psf_cf, int SigNx, Bool PsfMaxShift) +{ + FFTN_1D FFT1D; + int Nx = InPsf.nx(); + if (Nx < SigNx) Nx = SigNx; + + int Nx1=Nx; + // Bool ModifSize=False; + // dec_line_column (Nl, Nl1); + // dec_line_column (Nc, Nc1); + // if (Nl1 < Nc1) Nl1 = Nc1; + // else Nc1 = Nl1; + // if ((Nl != Nl1) || (Nc != Nc1)) ModifSize = True; + + fltarray OutPsf(Nx1, "Psf"); + if (PsfMaxShift == True) dec_center_psf (InPsf, OutPsf); + else psf_extend (InPsf, OutPsf); + // INFO(InPsf,"IN PSF"); + // INFO(OutPsf,"OUT PSF"); + + cfarray OutPsf_cf(Nx1, "Psf_cf"); + for (int i=0;i= N) Val = 2 * (N - 1) - ind; + else Val = ind; + } + if ((Val >= N) || (Val < 0)) Val = -1; + return (Val); +} + +/**************************************************************************/ +/********** ********************************************/ +/********** sig_extend ********************************************/ +/********** ********************************************/ +/**************************************************************************/ +/* extend signal with size of signal out (centered) */ +/**************************************************************************/ +void sig_extend (const fltarray& Signal, cfarray& Signal_Out) { + + int Nx0 = Signal.nx(); + int Nx1 = Signal_Out.nl(); + int i1,i0,Depi; + complex_f Zero = complex_f(0., 0.); + + Depi = (Nx1 - Nx0) / 2; + + for (i1 = 0; i1 < Nx1; i1++) { + + i0 = test_indice (i1 - Depi, Nx0); + if (i0 < 0) Signal_Out(i1) = Zero; + // if ((i0 < 0) || (j0 < 0)) Signal_Out(i1) = Zero; + else Signal_Out (i1) = complex_f(Signal(i0), 0.); + } +} + + + +/**************************************************************************/ +/********** ********************************************/ +/********** sig_extract ********************************************/ +/********** ********************************************/ +/**************************************************************************/ +/* extract cenrtered signal of size Nx1 in Out Signal */ +/**************************************************************************/ +void sig_extract (const cfarray &Signal, fltarray &Signal_Out, float Norm) +{ + int Nx0 = Signal_Out.nx(); + int Nx1 = Signal.nx(); + int i1,i0,Depi; + + Depi = (Nx1 - Nx0) / 2; + + for (i0 = 0; i0 < Nx0; i0++) { + + i1 = i0 + Depi; + Signal_Out (i0) = Signal (i1).real() / Norm; + } +} + + + +/**************************************************************************/ +/********** *******************************************/ +/********** dec_inverse *******************************************/ +/********** *******************************************/ +/**************************************************************************/ +/* compute fft(sig)*conj(fft(psf) */ +/**************************************************************************/ +void dec_inverse (fltarray& Signal, cfarray& Psf_cf, fltarray& Result, + float Eps=0.001, float Regul=0.); +void dec_inverse (fltarray& Signal, cfarray& Psf_cf, fltarray& Result, + float Eps, float Regul) { + FFTN_1D FFT1D; + int i; + int Nx = Signal.nx(); + int Nx1 = Psf_cf.nx(); + cfarray O_cf (Nx1, "Buffer conv"); + // float FluxPsf = (FluxNorm == False) ? 1.: Psf_cf(Nl1/2,Nc1/2).real(); + float FluxPsf = 1. / Psf_cf(Nx1/2).real(); + + if (Nx != Nx1) sig_extend (Signal, O_cf); + else { + for (i=0;i Eps) O_cf(i) /= Den; + else O_cf(i) = complex_f(0.); + } + FFT1D.fftn1d(O_cf,True); + // fft2d (O_cf, -1); + if (Nx != Nx1) sig_extract (O_cf, Result, FluxPsf); + else for (i = 0; i < Nx; i++) Result(i) = O_cf(i).real() / FluxPsf; +} + + + +/**************************************************************************/ +/********** ********************************************/ +/********** psf_convol ********************************************/ +/********** ********************************************/ +/**************************************************************************/ +/* compute convolution between signal and psf (in signal out) */ +/**************************************************************************/ +void psf_convol (fltarray& Signal, cfarray& Psf_cf, fltarray& Signal_out, + Bool FluxNorm) { + + FFTN_1D FFT1D; + int i,Nx = Signal.nx(); + int Nx1 = Psf_cf.nx(); + cfarray O_cf (Nx1, "Buffer conv"); + float FluxPsf = (FluxNorm == False) ? 1.: Psf_cf(Nx1/2).real(); + double Flux; + + if (Nx != Nx1) sig_extend (Signal, O_cf); + else { + for (i = 0; i < Nx; i++) O_cf(i) = complex_f(Signal(i),0); + } + // fft2d (O_cf); + FFT1D.fftn1d(O_cf); + O_cf *= Psf_cf; + // fft2d (O_cf, -1); + FFT1D.fftn1d(O_cf,True); + if (Nx != Nx1) { + + sig_extract (O_cf, Signal_out, FluxPsf); + Flux = Signal.total() / Signal_out.total(); + for (i = 0; i < Nx; i++) + Signal_out(i) = (float) (Signal_out(i) * Flux); + } + else + for (i = 0; i < Nx; i++) + Signal_out(i) = O_cf(i).real() / FluxPsf; + + // INFO(Imag,"IN ima"); + // INFO(Imag_out,"OUT ima"); +} +void psf_convol (fltarray& Signal, cfarray& Psf_cf, Bool FluxNorm) { + psf_convol (Signal, Psf_cf, Signal, FluxNorm); +} +void psf_convol (fltarray& Signal, fltarray& Psf, Bool FluxNorm) { + cfarray Psf_cf; + psf_get (Psf, Psf_cf, Signal.nx(), True); + psf_convol (Signal, Psf_cf, Signal, FluxNorm); +} +void psf_convol (fltarray& Signal, fltarray& Psf, fltarray& Signal_out, + Bool FluxNorm) { + cfarray Psf_cf; + psf_get (Psf, Psf_cf, Signal.nx(), True); + psf_convol (Signal, Psf_cf, Signal_out, FluxNorm); +} + +/**************************************************************************/ +/********** **************************************/ +/********** psf_convol_conj **************************************/ +/********** **************************************/ +/**************************************************************************/ +void psf_convol_conj (fltarray& Signal, cfarray& Psf_cf, fltarray& Signal_out) { + + FFTN_1D FFT1D; + int Nx = Signal.nx(); + int Nx1 = Psf_cf.nx(); + cfarray O_cf (Nx1, "Buffer conv"); + float FluxPsf = Psf_cf(Nx1/2).real(); + + if (Nx != Nx1) sig_extend (Signal, O_cf); + else { + for (int i=0; ithreshold(MR1D_Data); + + if (KillLastScale == True) + { + if (AdjointRec == False) + { + int s = MR1D_Data.nbr_band()-1; + for(int w=0; w MemModel(i)) + Gradient(i) -= RegulParam*log(Obj(i)/MemModel(i)); + } + else + { + for (i=0; i FLOAT_EPSILON) + Gradient(i) -= RegulParam*(1+log(Obj(i))); + } + break; + case DEC_MAP: + case DEC_MR_MAP: + for (i=0; i FLOAT_EPSILON) + Gradient(i) = Resi(i) / Signal_n(i); + else + Gradient(i) = FLOAT_EPSILON; + psf_convol_conj (Gradient, Psf_cf); + for (i=0; i FLOAT_EPSILON) + Gradient(i) = Resi(i) / Signal_n(i); + else Gradient(i) = 0.; + } + psf_convol_conj (Gradient, Psf_cf); + Gradient *= Obj; + // if ((RegulParam > 0) && + // RSig.obj_regul (Obj, Gradient, RegulParam*Noise_Sig); + break; + case DEC_MARKOV: + Gradient = Resi; + psf_convol_conj (Gradient, Psf_cf); + if (RegulParam > 0) RSig.obj_regul(Obj, Gradient, RegulParam*Noise_Sig); + break; + case DEC_TIKHONOV: + Gradient = Resi; + psf_convol_conj (Gradient, Psf_cf); + if (RegulParam > 0) RSig.obj_regul(Obj, Gradient, RegulParam); + break; + case DEC_MR_GRADIENT: + case DEC_GRADIENT: + Gradient = Resi; + psf_convol_conj (Gradient, Psf_cf); + // if (RegulParam > 0) RSig.obj_regul(Obj, Gradient, RegulParam*Noise_Sig); + break; + case DEC_MR_CLEAN: + case DEC_MR_CITTERT: + case DEC_CITTERT: + Gradient = Resi; + // if (RegulParam > 0) + // RSig.obj_regul (Obj, Gradient, RegulParam*Noise_Sig); + break; + + default: + cerr << "mr_deconv: Not implemented in this procedure ... "; + cerr << endl; + exit (0); + break; + } +} + + +/**************************************************************************/ +/********** *****************************************/ +/********** sig_gaussian *****************************************/ +/********** *****************************************/ +/**************************************************************************/ + +fltarray sig_gaussian (int Nx, float Fwhm, int Indi=-1); +fltarray sig_gaussian (int Nx, float Fwhm, int Indi) +{ + int Delta_i; + float sigma,sigma2,Energ = 0.; + register int i; + fltarray *Result = new fltarray(Nx, "gaussienne"); + int Ci=Indi; + + if (Ci < 0) Ci = Nx / 2; + sigma = 0.5 * Fwhm / sqrt (2. * log ((double) 2.)); + sigma2 = -2. * sigma * sigma; + + for (i=0; i FLOAT_EPSILON)) { + if ((It == 0) || (ABS(x - Mean) < Sm)) { + S0 ++; S1 += x; S2 += x*x; + } + } + } + if (S0 == 0) S0=1; + Mean = S1 / S0; + Sigma = S2/S0- Mean * Mean; + if ( Sigma >= 0.) Sigma = sqrt(Sigma); + else Sigma = 0.; + Sm = 3. * Sigma; + } + Moy= (float) Mean; + Sig =(float) Sigma; +} + + + + +/**************************************************************************/ +/********** *****************************************/ +/********** convol_gauss *****************************************/ +/********** *****************************************/ +/**************************************************************************/ + +void MR1D_Deconv::convol_gauss (fltarray& Data, float FWHM) +{ + fltarray Gauss; + Gauss = sig_gaussian (Data.nx(), FWHM); + norm_flux (Gauss); + psf_convol (Data, Gauss, True); +} + + +/**************************************************************************/ +/********** **************************************/ +/********** kill_last_scale **************************************/ +/********** **************************************/ +/**************************************************************************/ +void kill_last_scale (MR_1D & MR1D_Data) { + + int s = MR1D_Data.nbr_band()-1; + for (int w=0; w < MR1D_Data.size_scale_np(s); w++) + MR1D_Data(s,w) = 0; +} + + +/**************************************************************************/ +/********** ******************************************/ +/********** init_deconv ******************************************/ +/********** ******************************************/ +/**************************************************************************/ +void MR1D_Deconv::init_deconv(fltarray* FirstGuess, fltarray* ICF) +{ + float Flux; + Nx = Signal.nx(); + cout << Nx << endl; + Bool NonAddMethod = False; + int i; + + // set NonAddMethod to True for LUCY and MAP + if ( (DecMethod == DEC_LUCY) || (DecMethod == DEC_MR_LUCY) + || (DecMethod == DEC_MAP) || (DecMethod == DEC_MR_MAP)) + NonAddMethod = True; + + // Multiplyer image allocation MEM in memthod + if ((DecMethod == DEC_MEM) || (DecMethod == DEC_MEM_MODEL)) + Mult.alloc(Nx,"mult"); + // if ((DecMethod == DEC_LUCY) || + // ((OptimParam == True) && (StatNoise == NOISE_POISSON))) + // KeepImagn = True; + if (KeepImagn == True) Signal_n.alloc(Nx,"Signal_n"); + + // if an ICF is introduced, it can be either a Gaussian + // or an image given by the user + // the PSF must be convolved by the ICF + if (GaussConv == True) convol_gauss(Psf, Fwhm); + else if (ICF != NULL) { + Sig_ICF.alloc(Psf_cf.nx(), "ICF"); + dec_center_psf (*ICF, Sig_ICF); + UseICF = True; + FFT1D.convolve(Psf, Sig_ICF, Psf); + } + + // centering the PSF (normally the max of PSF is at the center of the image) + // the size of Psf_cf is a power of 2, which is not necessary the size of + // Psf, or the image size + psf_get (Psf, Psf_cf, Nx, PsfMaxShift); + + // allocate Obj if FirstGuess==NULL + if (FirstGuess != NULL) + TypeInit = DEC_INIT_GUESS; + else + Obj.alloc(Nx, "object"); + + // initialize Obj in function of TypeInit (with 0, flat, Signal or guess) + switch (TypeInit) + { + case DEC_INIT_ZERO: Obj.init(); break; + case DEC_INIT_FLAT: + Flux = Signal.total() / (float) (Nx); + Obj.init(Flux); + break; + case DEC_INIT_IMA: + Obj = Signal; + break; + case DEC_INIT_GUESS: + Obj = *FirstGuess; + break; + } + + // allocate resi + Resi.alloc(Nx,"resi"); + + /* Noise estimation in the Data */ + if (StatNoise == NOISE_GAUSSIAN) + { + if (Noise_Sig < FLOAT_EPSILON) + { + Noise_Sig = detect_noise_from_med (Signal); + if (Verbose == True) + cout << "Sigma Noise = " << Noise_Sig << endl; + } + } + if (Noise_Sig < FLOAT_EPSILON) Noise_Sig = 1.; + + if ((WaveFilterResi == True) || (KillLastScale == True)) + { + + WaveFilterResi = True; + // MR1D_Data.alloc(Nl, Nc, ModelData->nbr_scale(), + // ModelData->type_trans(), "MR1D_Data"); + MR1D_Data.alloc (Nx, ModelData->type_trans(), ModelData->nbr_scale(), + ModelData->filter_bank(), ModelData->TypeNorm); + + MR1D_Data.Border = Border; + ModelData->model(Signal, MR1D_Data); + + if (KillLastScale == True) + { + MR1D_Data.transform(Obj); + if (TypeInit == DEC_INIT_IMA) ModelData->threshold(MR1D_Data); + kill_last_scale (MR1D_Data); + MR1D_Data.recons(Obj); + } + } + + // Zero value in the first guess leads to problem with + // Lucy and MAP methods + if (NonAddMethod == True) + { + for (i=0; i < Nx; i++) + { + if (Obj(i) < FLOAT_EPSILON) Obj(i) = FLOAT_EPSILON; + if (Signal(i) < FLOAT_EPSILON) Signal(i) = FLOAT_EPSILON; + } + } + switch (DecMethod) + { + case DEC_INVERSE: + case DEC_MARKOV_LUCY: + case DEC_MARKOV: + RSig.GradOperType = OPER_MARKOV_2; + break; + case DEC_TIKHONOV: + RSig.GradOperType = OPER_LAPLACIAN; + break; + case DEC_MR_MEM: + case DEC_MR_MEM_NOISE: + case DEC_MEM: + case DEC_MEM_MODEL: + case DEC_MAP: + break; + case DEC_CITTERT: + case DEC_GRADIENT: + case DEC_LUCY: + case DEC_CLEAN: + case DEC_MR_CITTERT: + case DEC_MR_GRADIENT: + case DEC_MR_LUCY: + case DEC_MR_MAP: + if (RegulParam > 0) ApplySoftRegul = True; + if (Verbose == True) cout << "Wavelet Regularization " << RegulParam << endl; + break; + case DEC_MR_CLEAN: + case DEC_MR_VAGUELET: + case DEC_MR_INVERSE: + break; + } +} + + +/**************************************************************************/ +/********** ***********************************/ +/********** find_optim_poisson ***********************************/ +/********** ***********************************/ +/**************************************************************************/ +float MR1D_Deconv::find_optim_poisson (fltarray& Gradient) { + + float Cvgparam = 1., OldCvgparam; + fltarray Delta(Nx,"Buff"); + int i,Iter=0; + int Np = Nx; + float MinVal=1., MaxVal=100.; + float Val, Deriv1=0., Deriv2=0.; + if (RegulParam > 0.5) MinVal = 1. / (2.*RegulParam); + + // sear the maximum value + for (i=0; i < Np; i++) + if (Gradient(i) < - FLOAT_EPSILON) + { + Val = - Obj(i) / Gradient(i); + if (Val < MaxVal) MaxVal = Val; + } + // cout << "interval: Min = " << MinVal << " Max = " << MaxVal << endl; + + if (MaxVal > 0.) { + + // dec_convol (Gradient, Psf_cf, Delta); + psf_convol (Gradient, Psf_cf, Delta, True); + // Cvgparam = (MinVal + MaxVal) / 2.; + do { + + Iter++; + OldCvgparam = Cvgparam; + // Calculate first and second derivative + for (i=0; i FLOAT_EPSILON) + Deriv1 += Signal(i)*Delta(i) / Val ; + + // Second derivative + Val *= Val; + if (ABS(Val) > FLOAT_EPSILON) + Deriv2 -= Signal(i)*Delta(i)*Delta(i) / Val; + } + Cvgparam -= Deriv1/ Deriv2; + + // cout << "Iter " << Iter << " Cvgparam = " << Cvgparam << endl; + } while ((Iter < 10) && ((ABS(Cvgparam-OldCvgparam) > 0.1*OldCvgparam))); + if (Cvgparam > MaxVal) Cvgparam = MaxVal; + } + if (Cvgparam < MinVal) Cvgparam = MinVal; + return Cvgparam; +} + + + +/**************************************************************************/ +/********** **********************************/ +/********** find_optim_tikhonov **********************************/ +/********** **********************************/ +/**************************************************************************/ +float MR1D_Deconv::find_optim_tikhonov(fltarray& Gradient) { + + float Cvgparam = IterCvg; + fltarray Buff(Nx,"Buff"); + float LapG,LapO, Num, Den; + int i; + + if (RegulParam > 0.) { + + // dec_convol (Gradient, Psf_cf, Buff); + psf_convol(Gradient, Psf_cf, Buff, True); + Num = (Buff*Resi).total(); + Den = Buff.energy(); + for (i=0; i 0.5) MinVal = 1. / (2.*RegulParam); + // dec_convol (Gradient, Psf_cf, Buff); + psf_convol(Gradient, Psf_cf, Buff, True); + Cvgparam = (Buff*Resi).total()/Buff.energy(); + if (Cvgparam < MinVal) Cvgparam = MinVal; + else if (Cvgparam > 10) Cvgparam = 10.; + return Cvgparam; +} + + + +/**************************************************************************/ +/********** ****************************************/ +/********** im_find_optim ****************************************/ +/********** ****************************************/ +/**************************************************************************/ +float MR1D_Deconv::im_find_optim (fltarray& Gradient) { + + float Cvgparam = 1.; + if (DecMethod == DEC_TIKHONOV) Cvgparam = find_optim_tikhonov(Gradient); + else if (StatNoise == NOISE_GAUSSIAN) Cvgparam = find_optim_xi2(Gradient); + else if (StatNoise == NOISE_POISSON) Cvgparam = find_optim_poisson(Gradient); + + if (Verbose == True) { + cout << "Optim Val = " << Cvgparam << endl; + } + + if (Cvgparam < 0) Cvgparam = 0.; + return Cvgparam; +} + + + +/**************************************************************************/ +/********** *******************************************/ +/********** fonctional *******************************************/ +/********** *******************************************/ +/**************************************************************************/ +float MR1D_Deconv::fonctional() { + + float Val,Func = 0.; + int i; + + // Maxumum likehood part of the functional + if (StatNoise == NOISE_GAUSSIAN) Func = Resi.energy(); + else { // POISSON noise + + Func = Signal.energy(); // we add a constant in order to have + // a positive functional + for (i=0; i 0) + Func -= Signal(i)*log(Signal_n(i)); + } + + // regularized part of the functional + switch (DecMethod) + { + case DEC_TIKHONOV: + for (i=0; i 0) + Func += RegulParam*Obj(i)*log(Obj(i)); + break; + case DEC_MEM_MODEL: + if (UseModel == False) { + for (i=0; i 0) + Func += RegulParam*Obj(i)*log(Obj(i)); + } else { + for (i=0; i 0) + Func += RegulParam*((-Obj(i) + MemModel(i)) + + Obj(i)*log(Obj(i)/MemModel(i))); + } + break; + case DEC_MARKOV_LUCY: + case DEC_MARKOV: + case DEC_MR_MEM: + case DEC_MR_MEM_NOISE: + case DEC_MAP: + case DEC_CITTERT: + case DEC_GRADIENT: + case DEC_INVERSE: + case DEC_LUCY: + case DEC_CLEAN: + case DEC_MR_CITTERT: + case DEC_MR_GRADIENT: + case DEC_MR_LUCY: + case DEC_MR_MAP: + case DEC_MR_CLEAN: + case DEC_MR_VAGUELET: + case DEC_MR_INVERSE: + break; + } + + return Func; +} + + + +/**************************************************************************/ +/********** ***************************************/ +/********** im_iter_deconv ***************************************/ +/********** ***************************************/ +/**************************************************************************/ + +void MR1D_Deconv::sig_iter_deconv() +{ + fltarray Gradient(Nx,"Info"); + float Func, OldFunc=1e10, Sigma; + int i,Iter=0; + Bool Stop = False; + float Cvgparam = IterCvg; + float Sigma2N = Noise_Sig*Noise_Sig*Nx; + if (NormFlux == True) FluxImag = Signal.total(); + else FluxImag=0.; + + // Residual estimation + compute_resi(); + // INFO(Resi, "RESI"); + // start iterations + if ((Verbose == True) && (RegulParam > 0)) + cout << "Start iterate: Regul = " << RegulParam << endl; + do { + + // compute the gradient + sig_grad(Gradient); + + if (OptimParam == True) Cvgparam = im_find_optim(Gradient); + + // calculate the next object estimation with positivity constraint + if ((DecMethod == DEC_MAP) || (DecMethod == DEC_MR_MAP)) { + for (i = 0; i < Nx; i++) + Obj(i) *= Gradient(i); + } else if (DecMethod == DEC_MEM) + Obj = Gradient; + else { + for (i=0; i= 0) && (Iter % 1 == 0))) { + + Func = fonctional() / (float) (Nx); + printf("%d: Sigma = %f, Func = %f\n", Iter, Sigma, Func); + } + + Iter ++; + if (Iter > MaxIter) Stop = True; + else if (EpsCvg > FLOAT_EPSILON) { + + if (OldFunc < Sigma) Stop = True; + if ((ABS(OldFunc - Sigma)) < EpsCvg) Stop = True; + } + OldFunc = Sigma; + } while (Stop == False); + + if (Verbose == True) + { + Obj.info ("Solution: "); + Resi.info ("Resi: "); + } +} + +/**************************************************************************/ +/********** *******************************************/ +/********** sig_deconv *******************************************/ +/********** *******************************************/ +/**************************************************************************/ +void MR1D_Deconv::sig_deconv (fltarray* FirstGuess, fltarray* ICF) { + + + float Mean,Sigma; // used in DEC_CLEAN + float GammaClean = RegulParam; // used in DEC_CLEAN + Bool IterMethod = True; + + switch (DecMethod) + { + case DEC_INVERSE: + IterMethod = False; + StatNoise = NOISE_GAUSSIAN; + init_deconv(FirstGuess, ICF); + dec_inverse (Signal, Psf_cf, Obj, EpsCvg); + compute_resi(); + break; + case DEC_TIKHONOV: + case DEC_MARKOV: + case DEC_GRADIENT: + StatNoise = NOISE_GAUSSIAN; + break; + case DEC_MAP: + TypeInit = DEC_INIT_FLAT; + OptimParam = False; + StatNoise = NOISE_POISSON; + break; + case DEC_MARKOV_LUCY: + case DEC_LUCY: + TypeInit = DEC_INIT_FLAT; + StatNoise = NOISE_POISSON; + break; + case DEC_CITTERT: + case DEC_MEM: + case DEC_MEM_MODEL: + TypeInit = DEC_INIT_FLAT; + NormFlux = True; + OptimParam = False; + StatNoise = NOISE_GAUSSIAN; + break; + case DEC_CLEAN: + IterMethod = False; + StatNoise = NOISE_GAUSSIAN; + Nx = Signal.nx(); + Obj.alloc(Nx, "Obj"); + Resi.alloc(Nx, "Resi"); + Resi = Signal; + norm_flux (Psf, 1.); + if (Noise_Sig < FLOAT_EPSILON) + Noise_Sig = detect_noise_from_med (Resi); + sigma_clip (Resi,Mean,Sigma); + cout << "Mean+N_Sigma*Noise_Sig = " << Mean+N_Sigma*Noise_Sig << endl; + cout << "GammaClean = " << GammaClean << endl; + if (GaussConv == True) cout << "GaussConv = " << Fwhm << endl; + + dec_clean (Resi, Psf, Obj, 0., Mean+N_Sigma*Noise_Sig, + GammaClean, MaxIter, + DEFAULT_CLEAN_KAVER, DEFAULT_CLEAN_ADDRESI, + True, False, True); + if (ICF != NULL) { + Sig_ICF.alloc(Nx, "Psf"); + dec_center_psf (*ICF, Sig_ICF); + UseICF = True; + } + break; + case DEC_MR_VAGUELET: + IterMethod = False; + WaveFilterResi = True; + init_deconv(FirstGuess,ICF); + vaguelet(); + compute_resi(); + break; + case DEC_MR_CITTERT: + case DEC_MR_GRADIENT: + WaveFilterResi = True; + // NormFlux = True; + break; + case DEC_MR_LUCY: + case DEC_MR_MAP: + TypeInit = DEC_INIT_FLAT; + WaveFilterResi = True; + // NormFlux = True; + break; + case DEC_MR_CLEAN: + case DEC_MR_MEM: + case DEC_MR_MEM_NOISE: + case DEC_MR_INVERSE: + default: + cerr << "mr_deconv: Not implemented in this procedure ... "; + cerr << endl; + exit (0); + break; + } + if (IterMethod == True) + { + init_deconv(FirstGuess,ICF); + sig_iter_deconv(); + } + + if (GaussConv == True) convol_gauss(Obj, Fwhm); + else if (UseICF == True) psf_convol(Obj, Sig_ICF, Obj, True); + +} + +/************************************************************************/ + +void MR1D_Deconv::va_inverse_data(fltarray &NoiseSimu) +{ + if (StatNoise == NOISE_GAUSSIAN) simu_gaussian_noise(NoiseSimu, Noise_Sig); + else + { + // Find the noise map by filtering the data + MR1DFiltering CFilter(*ModelData, FIL_1D_THRESHOLD); + CFilter.Verbose = Verbose; + CFilter.filter(Signal, Obj); + + for (int i=0; i < Nx; i++) NoiseSimu(i) = Signal(i) - Obj(i); + } + + // inverse deconvolution + // dec_std_inverse (Imag, Psf_cf, Obj); + // dec_std_inverse (NoiseSimu, Psf_cf, NoiseSimu); + dec_inverse(Signal, Psf_cf, Obj, EpsCvg,RegulParam); + dec_inverse(NoiseSimu, Psf_cf, NoiseSimu, EpsCvg,RegulParam); +} + +/************************************************************************/ + +void MR1D_Deconv::vaguelet() +{ + int b,i; + int Nb = MR1D_Data.nbr_band()-1; + fltarray ThresholdMax(Nb+1); + fltarray NoiseSimu(Nx); + fltarray Scale; + + // inverse deconvolution on both imag and noise map + va_inverse_data(NoiseSimu); + // io_write_ima_float("xx_inv.fits", Obj); + // io_write_ima_float("xx_noise.fits", NoiseSimu); + + fltarray TabSigma(Nb); + MR1D_Data.transform(NoiseSimu); + for (b=0; b < Nb; b++) + { + MR1D_Data.scale (Scale, b); + ThresholdMax(b) = Scale.sigma() * ModelData->NSigma[b]; + } + for (b=0; b < Nb; b++) cout << "Band " << b << " sigma = " << TabSigma(b) << " T = " << ThresholdMax(b) << endl; + + // Update the multiresolution support + MR1D_Data.transform(Obj); + for (b=0; b < Nb; b++) + for (i=0; i < MR1D_Data.size_scale_np(b); i++) + { + if (ModelData->OnlyPositivDetect == True) + { + if (MR1D_Data(b,i) > ThresholdMax(b)) ModelData->support(b,i) = 1; + } + else + { + if (ABS(MR1D_Data(b,i)) > ThresholdMax(b)) ModelData->support(b,i) = 1; + } + } + + type_1d_filter Filter = FIL_1D_THRESHOLD; + MR1DFiltering CFilter(*ModelData, Filter); + NoiseSimu=Obj; + CFilter.Max_Iter = MaxIter; + CFilter.Epsilon = EpsCvg; + CFilter.Sup_Set = False; + CFilter.Verbose = Verbose; + CFilter.KillLastScale = KillLastScale; + CFilter.filter(NoiseSimu, Obj); // mr_support_filter(NoiseSimu, MR1D_Data, Obj); + if (PositivConstraint == True) + for (i=0; i