/*
                           COPYRIGHT (c) 1990
                     Kapteyn Astronomical Institute
                University of Groningen, The Netherlands
                           All Rights Reserved.

#>             editset.dc1

Program:       EDITSET

Purpose:       Program to edit data in a set.

Category:      MANIPULATION

File:          editset.c

Author:        M. Vogelaar

Keywords:

   INSET=      Give set (, subsets) to edit:

               Maximum number of subsets is 2048.
               You will be asked in a loop to enter a box for which
               you want to edit your data. The loop is aborted with
               carriage return at the ANOTHER= keyword.


   BOX=        Frame for input subsets.                 [entire subset]
   
               Select a box in which all data must be changed.
               If you want to change only one pixel, then enter this
               position as a box with sizes 1 in both directions. If 
               you want to perform the same actions on all data in 
               BOX= then enter a valid expression at EXPRESSION=. 
               If you want to edit the pixels manually or want to use 
               the Hermes database and file functions ('descr', 'table',
               'image' & 'file') then switch to 'manual mode' at the 
               EXPRESSION= keyword.
                                            

   EXPRESSION= Give expression:                           [Manual mode]
                  
               Enter an expression to change data in the given
               box. The parameters are the axis names like RA and DEC
               and the parameter 'DATA', which stands for the old.
               The syntax of the expressions is described below.
               An example is: EXPRESSION=DATA*ABS(RA) which will 
               multiply all values in BOX= with the absolute value
               of the grid in RA. 
               If carriage return is pressed at EXPRESSION=, then Manual
               Mode is selected and the keyword NEWVAL= is prompted for
               all pixels in BOX=.


   NEWVAL=     New value at: (x1,..,xn)                     [old value]
   
               No expression was entered at EXPRESSION= and you are 
               prompted with NEWVAL= and the position & value of the 
               current pixel starting at the lower left position in BOX=.
               It is possible to change the value of just one pixel or
               to change the values of a series of pixels at once
               (with a maximum of 250). This keyword recognizes the
               Hermes database and file functions. 
               Example: NEWVAL=image(AURORA FREQ 0,-5 0 5 0)
               which will replace all values in BOX= by values found
               in subset AURORA FREQ 0 in the box -5 0 5 0.


   ANOTHER=    Do you want another box?                           Y/[N]
   
               You can repeat the actions described above, for
               different boxes. Therefore a loop over the BOX= keyword
               is made. This loop is aborted by ANOTHER=N




Description:   The program EDITSET is an image editor. It can replace
               pixel values by user given values, or by new values 
               generated by a mathematical expression.
               Example: for a RA-DEC-FREQ set, this expression can contain 
               the parameters: RA, DEC, FREQ and DATA.
               RA is the RA grid position of the current pixel, DEC
               the DEC position etc. and DATA is the value of the
               current pixel.
               The first pixel to be changed is the first pixel in BOX=
               along the first axis in the subset (RA). The order
               of evaluation is: increase grid position along RA
               for all RA axes along DEC.
               If x and y are arbitrary parameters, the expression
               may contain:

              1) functions; syntax ff(..); where ff is one of the
                 following available functions:
                 abs(x)         absolute value of x
                 sqrt(x)        square root of x
                 sin(x)         sine of x
                 asin(x)        inverse sine of x
                 cos(x)         cosine of x
                 acos(x)        inverse cosine of x
                 tan(x)         tangent of x
                 atan(x)        inverse tan of x
                 atan2(x,y)     inverse tan (mod 2pi) x = sin, y = cos
                 exp(x)         exponential of x
                 ln(x)          natural log of x
                 log(x)         log (base 10) of x
                 sinh(x)        hyperbolic sine of x
                 cosh(x)        hyperbolic cosine of x
                 tanh(x)        hyperbolic tangent of x
                 rad(x)         convert x to radians
                 deg(x)         convert x to degrees
                 erf(x)         error function of x
                 erfc(x)        1-error function
                 max(x,y)       maximum of x and y
                 min(x,y)       minimum of x and y
                 sinc(x)        sin(x)/x (sinc function)
                 sign(x)        sign of x (-1,0,1)
                 mod(x,y)       gives remainder of x/y
                 int(x)         truncates to integer
                 nint(x)        nearest integer
                 ranu(x,y)      generates uniform noise between x and y
                 rang(x,y)      generates gaussian noise with mean x
                                and dispersion y
                 ranp(x)        generates poisson noise with mean x
                 ifeq(x,y,a,b)  returns a if x == y, else b
                 ifne(x,y,a,b)  returns a if x != y, else b
                 ifgt(x,y,a,b)  returns a if x > y, else b
                 ifge(x,y,a,b)  returns a if x >= y, else b
                 iflt(x,y,a,b)  returns a if x < y, else b
                 ifle(x,y,a,b)  returns a if x <= y, else b
                 ifblank(x,a,b) returns a if x == BLANK, else b
              2) constants; syntax cc; where cc is one of the following
                 available constants:
                 PI             3.14159....
                 C              speed of light (SI)
                 H              Planck (SI)
                 K              Boltzmann (SI)
                 G              gravitation (SI)
                 S              Stefan-Boltzman (SI)
                 M              mass of sun (SI)
                 P              parsec (SI)
                 BLANK          Universal undefined value
                 Note: the Hubble constant is not included.
              3) operators; syntax op; where op is one of the following
                 available operators:
                 +              addition
                 -              subtraction
                 *              multiplication
                 /              division
                 **             power



Example:       Multiply pixels in selected box by the absolute
               value of its RA grid position:

               EXPRESSION=DATA*ABS(RA)

Updates:       Mar 12, 1991: VOG, Document created.
               Feb 24, 1993, VOG, Changed default for BOX=


#<

*/

#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "math.h"
#include "cmain.h"
#include "gipsyc.h"
#include "init.h"
#include "finis.h"
#include "gdsinp.h"
#include "gdsc_ndims.h"
#include "setfblank.h"
#include "myname.h"
#include "anyout.h"
#include "nelc.h"
#include "cancel.h"
#include "gdsc_range.h"
#include "gdsc_grid.h"
#include "gdsbox.h"
#include "gdsc_fill.h"
#include "gdsi_read.h"
#include "gdsi_write.h"
#include "userint.h"
#include "usertext.h"
#include "userlog.h"
#include "fieini.h"
#include "fiedo.h"
#include "fiepar.h"
#include "cancel.h"
#include "gdsd_rchar.h"
#include "error.h"
#include "stabar.h"
#include "gdsc_name.h"
#include "userfio.h"
#include "userreal.h"
#include "userlog.h"
#include "userint.h"
#include "reject.h"


#define AXESMAX    10               /* Max. allowed number of axes in a set */
#define SUBSMAX    2048             /* Max. number of substructures to be specified */
#define MAXBUF     4096             /* Buffer size for I/O */
#define BIGSTORE   2048             /* Length of a string */
#define MAXBINS    500              /* Max. number of histogram bins */
#define VERSION    "1.0"            /* Version number of this program */
#define NONE       0                /* Default values for use in userxxx routines */
#define REQUEST    1
#define HIDDEN     2
#define EXACT      4
#define MAXREALS   250              /* Current maximum items for 'userreal' */
#define NO         0
#define YES        1


/* Keywords and messages */

#define KEY_INSET         tofchar("INSET=")
#define MES_INSET         tofchar("Give set (, subsets) to edit: " )
#define KEY_BOX           tofchar("BOX=")
#define MES_BOX           tofchar(" ")
#define KEY_NEWVAL        tofchar("NEWVAL=")  /* Has to generate message */
#define KEY_EXPRESSION    tofchar("EXPRESSION=")
#define MES_EXPRESSION    tofchar("Give expression:         [Manual mode]")
#define KEY_ANOTHER       tofchar("ANOTHER=")
#define MES_ANOTHER       tofchar("Do you want another box?    Y/[N]")


/* Initialize string with macro */
#define fmake(fchr,size) { \
                            static char buff[size+1]; \
                            int i; \
                            for (i = 0; i < size; buff[i++] = ' '); \
                            buff[i] = 0; \
                            fchr.a = buff; \
                            fchr.l = size; \
                         }


/* Malloc version of 'fmake'  */
#define finit( fc , len ) { fc.a = malloc( ( len + 1 ) * sizeof( char ) ) ;  \
                            fc.a[ len ] = '\0' ; \
                            fc.l = len ; }


#define MYMAX(a,b) ((a) > (b) ? (a) : (b))
#define MYMIN(a,b) ((a) > (b) ? (b) : (a))


/* Input of set, subsets: */

static fchar    Setin;                /* Name of the set */
static fint     subin[SUBSMAX];       /* Array for the subset coordinate words */
static fint     nsubsI;               /* Number of input subsets */
static fint     axnum[AXESMAX];       /* GDSINP axis numbers array */
static fint     axcount[AXESMAX];     /* GDSINP axis lengths array */
static fint     maxaxes  = AXESMAX;   /* Convert parameters to variables */
static fint     maxsubs  = SUBSMAX;
static int      subnr;                /* Index of current subset */
static int      i, m;                 /* Counters */
static fint     setlevel = 0;         /* Indicate set level */


/* Input of area etc.:*/

static fint     cwlo, cwhi;           /* Coordinate words */
static fint     FgridLO[AXESMAX];     /* Coordinate words for frame */
static fint     FgridHI[AXESMAX];
static fint     BgridLO[AXESMAX];     /* Coordinate words for box */
static fint     BgridHI[AXESMAX];


/* Data transfer: */

static fint     totpixels;            /* Total number of pixels in input */
static float    imageIN[MAXBUF];      /* Contains data to be changed */


/* Stabar related: */

static float    SBbegin;
static float    SBend;
static float    SBcurrent;


/* Expression evaluation: */


static float    parmvalues[MAXBUF*(AXESMAX+1)];
                                       /* There can be as much as          */
                                       /* AXESMAX+1 parameters in an       */
                                       /* expression. Those are the axis   */
                                       /* names and the parameter 'DATA'.  */
                                       /* For each parameter there is      */
                                       /* MAXBUF variable storage.         */


/* Miscellaneous: */

static float    blank;                 /* Value of system blank */
static char     messbuf[BIGSTORE];     /* Buffer for text message */
static int      substruct[AXESMAX+1];  /* Number of pixels in substructure with dimension n */
static char     posstr[BIGSTORE];      /* Contains converted position */
static int      pixelnr;               /* pixel pos. in 1-dim array */
static fint     gridpos[AXESMAX+1];
static bool     another;



static void numtoposstr( fchar Setin,
                         int   pixelnr,
                         int   *substruct,
                         fint  *BgridLO,
                         int   subdim,
                         char  *posstr )
/*------------------------------------------------------------*/
/* Convert an one dimensional pixel position in a box with    */
/* dimension n, to a n-dimensional pixel position. Put results*/
/*  in a string.                                              */
/*------------------------------------------------------------*/
{
   int     position[AXESMAX+1];
   int     i;
   char    dumbuf[20];
   fint    err;
   fint    grid;
   fint    setdim;


   setdim  = gdsc_ndims_c( Setin, &setlevel );
   strcpy( posstr, "(" );
   if (subdim == 0)
   {
      /* In the set specification, only a pixel is defined */
      for (i = 0; i < (int) setdim; i++)
      {
         /* Items Setin, axnum, and subin are global variables */
         err = 0;
         grid = gdsc_grid_c( Setin, &axnum[i], subin, &err );
         sprintf( dumbuf, "%-d", (int) grid );
         strcat( posstr, dumbuf );
         if (i < (int) setdim-1)
            strcat( posstr, "," );
      }
   }
   else
   {
      for (i = subdim; i >= 1; i--)
      {
         position[i] = pixelnr / ( substruct[i-1] );
         position[i] += (int) BgridLO[i-1];             /* Offset from corner */
         /* subtract number of pixels in underlying substructures */
         pixelnr -= (pixelnr / substruct[i-1] ) * substruct[i-1];
      }
      for (i = 1; i <= subdim; i++)
      {
         sprintf( dumbuf, "%-d", position[i] );
         strcat( posstr, dumbuf );
         if (i < subdim)
            strcat( posstr, "," );
      }
   }
   strcat( posstr, ")" );
}



static void numtopos( int    pixelnr,
                      int    *substruct,
                      fint   *BgridLO,
                      int    setdim,
                      int    subdim,
                      fint   *gridpos )
/*------------------------------------------------------------*/
/* Convert an one dimensional pixel position in a box with    */
/* dimension n, to a n-dimensional pixel position. Put results*/
/* in an array. The grid that belongs to the first axis, is   */
/* the first element in that array, but its index is 1!       */
/*------------------------------------------------------------*/
{
   int  i;
   fint err;


   if (subdim == 0)
   {
      /* In the set specification, only a pixel is defined */
      for (i = 0; i < (int) subdim; i++)
      {
         err = 0;
         /* Items Setin, axnum, and subin are global variables */
         gridpos[i] = gdsc_grid_c( Setin, &axnum[i], subin, &err );
      }
   }
   else
   {
      /* These axes are the repeat axes, specified by the user */
      for (i = setdim; i > subdim; i--)
      {
         err = 0;
         gridpos[i] = gdsc_grid_c( Setin, &axnum[i-1], subin, &err );
      }
      for (i = subdim; i >= 1; i--)
      {
         gridpos[i] = (fint) ( pixelnr / ( substruct[i-1] ) );
         gridpos[i] += (fint) BgridLO[i-1];         /* Offset from corner */
         /* subtract number of pixels in underlying substructures */
         pixelnr -= ( pixelnr / substruct[i-1] ) * substruct[i-1];
      }
   }
}




fint askexpression( fint   *Ffieid,
                    fchar  Setin,
                    fint   *axnum )
/*------------------------------------------------------------*/
/* Ask user for an valid expression. Determine all allowed    */
/* parameters (axis names and "DATA") and check the           */
/* expression. Return the number of parameters in the         */
/* expression.                                                */
/*------------------------------------------------------------*/
{
   fint   Fres;
   int    agreed;
   fchar  Fiestr;                                /* User specified expression */
   fint   errpos;                     /* Error position in a wrong expression */
   char   messbuf[BIGSTORE+26];
   fint   dfault;
   fint   numpars = -1;                   /* Number of parameters in an expr. */
   fchar  parnames;                           /* All possible parameter names */
   fint   npars;                           /* Max. number of par. in an expr. */
   fint   err;
   fchar  Fctype;
   char   dummystr[20];
   int    n;
   fint   setlevel = 0;
   fint   setdim;


   setdim  = gdsc_ndims_c( Setin, &setlevel );
   fmake( parnames, 10*(AXESMAX+1) );
   fmake( Fctype, 10 );
   fmake( Fiestr, BIGSTORE );
   do
   {
      agreed = NO;
      dfault = REQUEST;
      Fres = usertext_c( Fiestr, &dfault, KEY_EXPRESSION, MES_EXPRESSION );
      if (Fres == 0)
         return( -1 );                          /* Stop expression activities */

      /* Determine all possible 'fie' variable names */
      strcpy( parnames.a, "data         " );          /* Spaces are important */
      for (n = 0; n < setdim; n++)
      {
         err = 0;
         gdsc_name_c( Fctype, Setin, &axnum[n], &err );
         /* Copy until space or hyphen! */
         sprintf( dummystr, "%-10s", strtok( Fctype.a, " -" ) );
         strncat( parnames.a, dummystr, 10 );
      }

      npars = setdim + 1;
      /* Tell fie which parameters can be expected */
      parnames.l = 10;
      Fres = fiepar_c( parnames, &npars );

      /* Initialize the expression */
      Fres = fieini_c( Fiestr, Ffieid, &errpos );
      if (Fres >= 0)
      {
         numpars = Fres;
         agreed = (numpars <= (setdim+1));
         if (!agreed)
         {
            anyoutf( 1, "Number of parameters greater than number of axes." );
         }
      }
      if (Fres == -1)
      {
         agreed = NO;
         sprintf( messbuf,
                 "Syntax error in expression at position %-d",
                  errpos );
         anyoutf( 1, messbuf );
      }
      if (Fres == -2)
      {
         agreed = NO;
         anyoutf( 1, "No storage space left." );
      }
      if (!agreed)
      {
         reject_c( KEY_EXPRESSION, tofchar("Invalid expression") );
      }
   }
   while (!agreed);

   sprintf( messbuf,
           "Accepted expression: %.*s",
            nelc_c( Fiestr ),
            Fiestr.a );
   anyoutf( 1, messbuf );
   return( numpars );
}



MAIN_PROGRAM_ENTRY
/*-------------------------------------------------------------*/
/* EDITSET main. The macro MAIN_PROGRAM_ENTRY replaces the     */
/* C-call main() to start the main body of your GIPSY          */
/* application. Variables defined as 'fchar' start with a      */
/* capital.                                                    */
/*-------------------------------------------------------------*/
{

   fint    dfault, nitems;
   fint    subdim, setdim;
   fint    scrnum;
   fint    class;
   fint    r1, r2;
   int     validexpression;


   init_c();                               /* contact Hermes */
   /* Task identification */
   {
      fchar   Task;                        /* Name of current task */
      fmake( Task, 20 );                   /* Macro 'fmake' must be available */
      myname_c( Task );                    /* Get task name */
      Task.a[nelc_c(Task)] = '\0';        /* Terminate task name with null char. */
      IDENTIFICATION( Task.a, VERSION );   /* Show task and version */
   }

   setfblank_c( &blank );
   fmake(Setin, 80);


   dfault = NONE;
   subdim = 0;
   scrnum = 3;
   class  = 1;
   nsubsI = gdsinp_c( Setin,
                      subin,
                      &maxsubs,
                      &dfault,
                      KEY_INSET,
                      MES_INSET,
                      &scrnum,
                      axnum,
                      axcount,
                      &maxaxes,
                      &class,
                      &subdim );
   setdim  = gdsc_ndims_c( Setin, &setlevel );
   /*--------------------------------------------------*/
   /* Determine the edges of this its frame            */
   /* (FgridLO/HI)                                     */
   /*--------------------------------------------------*/
   r1 = 0;
   gdsc_range_c( Setin, &setlevel, &cwlo, &cwhi, &r1 );
   r1 = r2 = 0;
   for (m = 0; m < (int) setdim; m++) {
      r1 = r2 = 0;
      FgridLO[m] = gdsc_grid_c( Setin, &axnum[m], &cwlo, &r1 );
      FgridHI[m] = gdsc_grid_c( Setin, &axnum[m], &cwhi, &r2 );
   }

   do
   /*--------------------------------------------------*/
   /* Prepare a box for INSET.                         */
   /*--------------------------------------------------*/
   {
      fint     boxopt = 0;
      fint     fie_id;

      dfault = REQUEST;
      scrnum = 3;
      gdsbox_c( BgridLO, BgridHI,
                Setin, subin,
                &dfault,
                KEY_BOX, MES_BOX,
                &scrnum,
                &boxopt );
      /* Count number of pixels in this substructure */
      totpixels = 1;
      for(m = 0; m < (int) subdim; m++)                     /* For one subset */
         totpixels *= (BgridHI[m] - BgridLO[m] + 1);
      totpixels *= nsubsI;                                 /* For all subsets */

      /*--------------------------------------------------*/
      /* Determine number of pixels of substructures with */
      /* lower dimensions than the subset dimension.      */
      /* These numbers are used in the 'numtopos' routines*/
      /* to determine the 'setdim' dimensional grid       */
      /* position of a pixel.                             */
      /*--------------------------------------------------*/
      substruct[0] = 1;                                          /* One pixel */
      for (i = 1; i <= (int) subdim; i++ )
         substruct[i] = (int) ( BgridHI[i-1]-BgridLO[i-1]+1 ) * substruct[i-1];


      /*--------------------------------------------------*/
      /* Ask to specify an expression. If <return> is     */
      /* pressed, the variable askexpression' will be set */
      /* to a negative value, else the number of specified*/
      /* parameters will be returned.                     */
      /*--------------------------------------------------*/

      validexpression = (askexpression(&fie_id, Setin, axnum) >= 0);


      if (validexpression)
      {
         /* Initialize 'stabar' function */
         SBbegin   = 0.0;
         SBend     = (float) totpixels;
         SBcurrent = 0.0;
         stabar_c( &SBbegin, &SBend, &SBcurrent );
      }

      for(subnr = 0; subnr < nsubsI; subnr++)
      /*--------------------------------------------------*/
      /* Loop over all specified subsets.                 */
      /*--------------------------------------------------*/
      {
         fint   tidIN, tidOUT;                          /* File transfer id's */
         fint   maxIObuf = MAXBUF;
         fint   pixelsdone;

         /* Make coordinate words for these corners */
         cwlo    = gdsc_fill_c( Setin, &subin[subnr], BgridLO );
         cwhi    = gdsc_fill_c( Setin, &subin[subnr], BgridHI );
         tidIN   = 0;
         tidOUT  = 0;
         pixelnr = 0;
         do
         {
            gdsi_read_c( Setin,
                         &cwlo, &cwhi,
                         imageIN,
                         &maxIObuf,
                         &pixelsdone,
                         &tidIN );
            /*--------------------------------------------------*/
            /* There is no expression ==> change manually       */
            /*--------------------------------------------------*/
            if (!validexpression)
            {
               i = 0;
               do
               {
                  /* Loop over all pixels in buffer */
                  nitems = MYMIN(pixelsdone, MAXREALS);
                  dfault = REQUEST;
                  numtoposstr( Setin, pixelnr, substruct,
                               BgridLO, subdim, posstr );
                  if (imageIN[i] == blank)
                     sprintf( messbuf, "New val. at %s:      [blank]", posstr );
                  else
                     sprintf( messbuf, "New val. at %s:      [%.4g]", posstr, imageIN[i] );
                  r1 = userreal_c( &imageIN[i],
                                   &nitems,
                                   &dfault,
                                   KEY_NEWVAL,
                                   tofchar(messbuf) );
                  cancel_c( KEY_NEWVAL );
                  if (r1 == 0)
                  {
                     i++;
                     pixelnr++;
                  }
                  else
                  {
                     i += (int) r1;
                     pixelnr += (int) r1;
                  }
               }
               while (i < (int) pixelsdone);
            }
            else
            /*--------------------------------------------------*/
            /* There is an expression ==> change by expression. */
            /* Fill the first 'pixelsdone' number of old values */
            /* (corresponding to the DATA parameter) in the     */
            /* first part of the array. The array is divided in */
            /* setdim+1 parts of length 'pixelsdone'. The other */
            /* parts are filled with the grid positions of the  */
            /* pixels.                                          */
            /*--------------------------------------------------*/
            {
               int j,m;
               for (j = 0; j < (int) pixelsdone; j++)
               {
                  parmvalues[j] = imageIN[j];
                  numtopos( pixelnr++, substruct, BgridLO,
                            setdim, subdim, gridpos );
                  for (m = 1; m <= (int) setdim; m++) {
                     parmvalues[j + m*pixelsdone] = gridpos[m];
                  }
               }
               r1 = fiedo_c( parmvalues, &pixelsdone, imageIN, &fie_id );
            }
            gdsi_write_c( Setin,
                          &cwlo, &cwhi,
                          imageIN,
                          &maxIObuf,
                          &pixelsdone,
                          &tidOUT );
            SBcurrent += (float) pixelsdone;
            if (validexpression)
               stabar_c( &SBbegin, &SBend, &SBcurrent );
         }
         while (tidIN != 0);                          /* All pixels examined? */
      }                                                  /* All subsets done? */
      cancel_c( KEY_BOX );
      cancel_c( KEY_EXPRESSION );
      another = toflog( 0 );
      dfault  = REQUEST;
      nitems  = 1;
      r1      = userlog_c( &another,
                           &nitems,
                           &dfault,
                           KEY_ANOTHER,
                           MES_ANOTHER );
      another = tobool( another );
      cancel_c( KEY_ANOTHER );
   }
   while( another );                                          /* Another box? */

   finis_c();                                                  /* Quit Hermes */
   return( EXIT_SUCCESS );                                    /* Dummy return */
}

