/*  cfill.c
   
    For given field Q with defined and undefined regions and 
    a target grid point (X,Y), return new field Q' that keeps Q
    defined ONLY a contiguous defined region containing the target,
    and sets all other regions to undefined. This is achieved
    through a simple fill operation, which fills the area within the 
    contour that defines the target contiguous region 
    starting from the grid point at (X,Y).
   
   *  Arguments expected:
   *    1) field Q with defined and undefined regions [expr, grid]
   *    2) target X within a defined region [positive int]
   *    3) target Y within a defined region [positive int]
   * 
   *  Returns:
   *    1) field Q but undef everywhere EXCEPT in the
   *       contiguous, defined region that contains the target


    Example use in a GrADS script:
        
        ;* define 500-hPa geopotential height
        'set lev 500'
        'define gph=geopt/9.8'

        ;* get minimum geopt location where a strong 
        ;* cyclone of interest is 
        'd aminlocx(gph,lon=0,lon=140,lat=55,lat=85)'
        x=subwrd(result,4);
        'd aminlocy(gph,lon=0,lon=140,lat=55,lat=85)'
        y=subwrd(result,4);

        ;* problem:   the cyclone of interest is in a domain
        ;*            where other cyclonic features exist.
        ;*            Can we isolate the cyclone from the
        ;*            other features?

        ;* solution:  use cfill() to define the height field 
        ;*            of the strong cyclone centered at 
        ;*            (x,y) within the cyclone's outer-most
        ;*            closed height contour at 544-dam
        
        'define cyclgph=cfill(maskout(gph,5440-gph),'x','y')'

        [...]


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

#include "grads.h"    

gaint cfill (struct gafunc *pfc, struct gastat *pst, struct gaudpinfo *pinfo) {
  gaint i,j,rc;
  struct gagrid *pgr;
  char *valu;
  int xo,yo,inc;
  const int MAXN=1000*1000;
  gaint inew[MAXN];
  gaint jnew[MAXN];
  int Nnew;
  int io,jo,ID;
  char msg[126];

  size_t sz;
  struct gastat *tarst; // target field structure
  struct gagrid *targr; // target grid structure
  char *taru;           // target array def mask (0=undef)

  /* Check if this plug-in is still compatible with the GrADS interface */
  if (pinfo->version != UDPVERS) printf("Error: version mismatch. Recompile this plug-in\n");

  /* Check for the proper number of arguments */
  if (pfc->argnum!=3) {
    gaprnt (0,"Error from CFILL: Three arguments expected. \n");
    return (1);
  }

  sz = sizeof(struct gastat);
  tarst = (struct gastat *)galloc(sz,"funcpst");
  *tarst = *pst;
  
  /* Invoke gaexpr to evaluate the first argument, 
   * which should be a GrADS expression.
     The result gets stored in the gastat structure 'pst' */
  rc = gaexpr(pfc->argpnt[0],pst);
  if (rc) return (rc);

  /* evaluate 1st argument again (field to regrid) */
  rc = gaexpr(pfc->argpnt[0],tarst);
  if (rc) return (rc);
  
  /* If the expression is a grid ... */
  if (pst->type==1) {
    /* Set up pointers to the gridded data structure, the grid of data, and the undef mask  */
    pgr   = pst->result.pgr;
    valu  = pgr->umask;
    targr = tarst->result.pgr;
    taru  = targr->umask;

    /* let's convert arguments 2, 3
     * from strings to integers. */
    xo=0;
    i=0; while(pfc->argpnt[1][i] != '\0') {
      switch(pfc->argpnt[1][i]) {
        case '0': inc = 0;  break;
        case '1': inc = 1;  break;
        case '2': inc = 2;  break;
        case '3': inc = 3;  break;
        case '4': inc = 4;  break;
        case '5': inc = 5;  break;
        case '6': inc = 6;  break;
        case '7': inc = 7;  break;
        case '8': inc = 8;  break;
        case '9': inc = 9;  break;
        default:  inc = -1; break;
      }
      if (inc < 0) break;
      xo = (xo * 10) + inc;
      i++;
    }
    if (xo == 0) {
      gaprnt (0,"Error from CFILL: xo (arg 2) must be positive. \n");
      return (1);
    }
    
    yo=0;
    i=0; while(pfc->argpnt[2][i] != '\0') {
      switch(pfc->argpnt[2][i]) {
        case '0': inc = 0;  break;
        case '1': inc = 1;  break;
        case '2': inc = 2;  break;
        case '3': inc = 3;  break;
        case '4': inc = 4;  break;
        case '5': inc = 5;  break;
        case '6': inc = 6;  break;
        case '7': inc = 7;  break;
        case '8': inc = 8;  break;
        case '9': inc = 9;  break;
        default:  inc = -1; break;
      }
      if (inc < 0) break;
      yo = (yo * 10) + inc;
      i++;
    }
    if (yo == 0) {
      gaprnt (0,"Error from CFILL: yo (arg 3) must be positive. \n");
      return (1);
    }

    /* initialize the result field's undef mask (all undef)  */
    for (j=0; j < pgr->jsiz ; j++) {
      for (i=0; i < pgr->isiz ; i++) {
        *(valu + i + (j * pgr->isiz)) = 0;
      }
    }

    /* convert the initial coordinates to result grid-relative    */
    /* EX: Your domain in GrADS is set to                         */
    /*     X: 201 to 301; Y: 1 to 101.                            */
    /*     A point within the masked contour you want to isolate  */
    /*     is at X = 250, Y = 51.                                 */
    /*     The result grid in this plug-in has dimensions of      */
    /*     isiz=101, jsiz=101. We offset the target point at      */
    /*     (X,Y) by the dimension mins to get the location of the */
    /*     target point relative to (i=0,j=0) in the result array */
    io = xo - pgr->dimmin[0];
    jo = yo - pgr->dimmin[1];

    /* try to initialize fill at (xo,yo) */
    Nnew=0;
    if (io >= 0 && jo >= 0 && io < pgr->isiz && jo < pgr->jsiz) {
      ID = io + (jo * pgr->isiz);
      if (*(taru + ID)) {
        *(valu + ID) = 1;
        inew[Nnew]=io; 
        jnew[Nnew]=jo;
        Nnew++;
      }
    }

    /* do fill */
    while (Nnew>0 && Nnew < MAXN - 4) {
      
      /* i,j of a newly added fill point */
      io = inew[Nnew-1];
      jo = jnew[Nnew-1];

      /* check the four points adjacent to this new point 
       * to see if they should be filled too */
      Nnew--; // as we're checking this new point now,
              // subtract it from the N of new points.
      
      /* left */
      i=io-1;
      j=jo;
      if (i >= 0) {
        ID = i + (j * pgr->isiz);
        // fill this point if it's defined and needs filling
        if (!*(valu + ID) && *(taru + ID)) {
          *(valu + ID) = 1; // fill
          inew[Nnew]=i; // add this point at (i,j)
          jnew[Nnew]=j; // to the list of new fill points
          Nnew++; // add to the N of new points.
        }
      }

      /* right */
      i=io+1;
      j=jo;
      if (i < pgr->isiz) {
        ID = i + (j * pgr->isiz);
        // fill this point if it's defined
        if (!*(valu + ID) && *(taru + ID)) {
          *(valu + ID) = 1; // fill
          inew[Nnew]=i; // add this point at (i,j)
          jnew[Nnew]=j; // to the list of new fill points
          Nnew++;
        }
      }

      /* up */
      i=io;
      j=jo+1;
      if (j < pgr->jsiz) {
        ID = i + (j * pgr->isiz);
        // fill this point if it's defined
        if (!*(valu + ID) && *(taru + ID)) {
          *(valu + ID) = 1; // fill
          inew[Nnew]=i; // add this point at (i,j)
          jnew[Nnew]=j; // to the list of new fill points
          Nnew++;
        }
      }
      
      /* down */
      i=io;
      j=jo-1;
      if (j >= 0) {
        ID = i + (j * pgr->isiz);
        // fill this point if it's defined
        if (!*(valu + ID) && *(taru + ID)) {
          *(valu + ID) = 1; // fill
          inew[Nnew]=i; // add this point at (i,j)
          jnew[Nnew]=j; // to the list of new fill points
          Nnew++;
        }
      }
    }

    if (Nnew>0) {
      gaprnt (2,"Error from CFILL: Fill perimeter too large!\n");
      gaprnt (2,"   Your Perimeter size >= ");
      snprintf(msg,125,"%d",Nnew);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      gaprnt (2,"   Max perimeter size  == ");
      snprintf(msg,125,"%d",MAXN);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      /* debug statements ...
      snprintf(msg,125,"%d",pgr->jsiz);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%g",pgr->rmin);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%g",pgr->rmax);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%g",*(pgr->ivals));
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%d",pgr->dimmin[0]);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%d",pgr->dimmin[1]);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%d",pgr->dimmax[0]);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      snprintf(msg,125,"%d",pgr->dimmax[1]);
      gaprnt (2,msg);
      gaprnt (2,"\n");
      */
    }
  } 
  
  /* ...the expression is station data */
  else {              
    gaprnt (0,"Error from CFILL: bad data type (non-grid). \n");
    return (1);
  }
  
  gree(tarst,"f404");

  /* finish without errors */
  return (0);
}