#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "fitsio2.h"
#include "region.h"

static int Pt_in_Poly( double x, double y, int nPts, double *Pts );

/*---------------------------------------------------------------------------*/
int ffrrgn( const char *filename,
            WCSdata    *wcs,
            SAORegion  **Rgn,
            int        *status )
/*  Read regions from a SAO-style region file and return the information     */
/*  in the "SAORegion" structure.  If it is nonNULL, use wcs to convert the  */
/*  region coordinates to pixels.  Return an error if region is in degrees   */
/*  but no WCS data is provided.                                             */
/*---------------------------------------------------------------------------*/
{
   char     *currLine;
   char     *namePtr, *paramPtr, *currLoc;
   char     *pX, *pY, *endp;
   long     allocLen, lineLen, hh, mm, dd;
   double   *coords = 0, X, Y, R, x, y, ss, xsave= 0., ysave= 0.;
   int      nParams, nCoords, negdec;
   int      i, done;
   FILE     *rgnFile;
   coordFmt cFmt;
   SAORegion *aRgn;
   RgnShape *newShape, *tmpShape;

   if( *status ) return( *status );

   aRgn = (SAORegion *)malloc( sizeof(SAORegion) );
   if( ! aRgn ) {
      ffpmsg("Couldn't allocate memory to hold Region file contents.");
      return(*status = MEMORY_ALLOCATION );
   }
   aRgn->nShapes    =    0;
   aRgn->Shapes     = NULL;
   if( wcs && wcs->exists )
      aRgn->wcs = *wcs;
   else
      aRgn->wcs.exists = 0;

   cFmt = pixel_fmt; /* set default format */

   /*  Allocate Line Buffer  */

   allocLen = 512;
   currLine = (char *)malloc( allocLen * sizeof(char) );
   if( !currLine ) {
      free( aRgn );
      ffpmsg("Couldn't allocate memory to hold Region file contents.");
      return(*status = MEMORY_ALLOCATION );
   }

   /*  Open Region File  */

   if( (rgnFile = fopen( filename, "r" ))==NULL ) {
      sprintf(currLine,"Could not open Region file %s.",filename);
      ffpmsg( currLine );
      free( currLine );
      free( aRgn );
      return( *status = FILE_NOT_OPENED );
   }
   
   /*  Read in file, line by line  */

   while( fgets(currLine,allocLen,rgnFile) != NULL ) {

      /*  Make sure we have a full line of text  */

      lineLen = strlen(currLine);
      while( lineLen==allocLen-1 && currLine[lineLen-1]!='\n' ) {
         currLoc = (char *)realloc( currLine, 2 * allocLen * sizeof(char) );
         if( !currLoc ) {
            ffpmsg("Couldn't allocate memory to hold Region file contents.");
            *status = MEMORY_ALLOCATION;
            goto error;
         } else {
            currLine = currLoc;
         }
         fgets( currLine+lineLen, allocLen+1, rgnFile );
         allocLen += allocLen;
         lineLen  += strlen(currLine+lineLen);
      }

      currLoc = currLine;
      if( *currLoc == '#' ) {

         /*  Look to see if it is followed by a format statement...  */
         /*  if not skip line                                        */

         currLoc++;
         while( *currLoc==' ' ) currLoc++;
         if( !strncasecmp( currLoc, "format:", 7 ) ) {
            if( aRgn->nShapes ) {
               ffpmsg("Format code encountered after reading 1 or more shapes.");
               *status = PARSE_SYNTAX_ERR;
               goto error;
            }
            currLoc += 7;
            while( *currLoc==' ' ) currLoc++;
            if( !strncasecmp( currLoc, "pixel", 5 ) ) {
               cFmt = pixel_fmt;
            } else if( !strncasecmp( currLoc, "degree", 6 ) ) {
               cFmt = degree_fmt;
            } else if( !strncasecmp( currLoc, "hhmmss", 6 ) ) {
               cFmt = hhmmss_fmt;
            } else if( !strncasecmp( currLoc, "hms", 3 ) ) {
               cFmt = hhmmss_fmt;
            } else {
               ffpmsg("Unknown format code encountered in region file.");
               *status = PARSE_SYNTAX_ERR;
               goto error;
            }
         }

      } else if( !strncasecmp( currLoc, "glob", 4 ) ) {
		  /* skip lines that begin with the word 'global' */

	  } else {

         while( *currLoc != '\0' ) {

            namePtr  = currLoc;
            paramPtr = NULL;
            nParams  = 1;

            /*  Search for closing parenthesis  */

            done = 0;
            while( !done && !*status && *currLoc ) {
               switch (*currLoc) {
               case '(':
                  *currLoc = '\0';
                  currLoc++;
                  if( paramPtr )   /* Can't have two '(' in a region! */
                     *status = 1;
                  else
                     paramPtr = currLoc;
                  break;
               case ')':
                  *currLoc = '\0';
                  currLoc++;
                  if( !paramPtr )  /* Can't have a ')' without a '(' first */
                     *status = 1;
                  else
                     done = 1;
                  break;
               case '#':
               case '\n':
                  *currLoc = '\0';
                  if( !paramPtr )  /* Allow for a blank line */
                     done = 1;
                  break;
               case ':':  
                  currLoc++;
                  cFmt = hhmmss_fmt;
                  break;
               case 'd':
                  currLoc++;
                  cFmt = degree_fmt;
                  break;
               case ',':
                  nParams++;  /* Fall through to default */
               default:
                  currLoc++;
                  break;
               }
            }
            if( *status || !done ) {
               ffpmsg( "Error reading Region file" );
               *status = PARSE_SYNTAX_ERR;
               goto error;
            }

            /*  Skip white space in region name  */

            while( *namePtr==' ' ) namePtr++;

            /*  Was this a blank line? Or the end of the current one  */

            if( ! *namePtr && ! paramPtr ) continue;

            /**************************************************/
            /*  We've apparently found a region... Set it up  */
            /**************************************************/

            if( !(aRgn->nShapes % 10) ) {
               if( aRgn->Shapes )
                  tmpShape = (RgnShape *)realloc( aRgn->Shapes,
                                                  (10+aRgn->nShapes)
                                                  * sizeof(RgnShape) );
               else
                  tmpShape = (RgnShape *) malloc( 10 * sizeof(RgnShape) );
               if( tmpShape ) {
                  aRgn->Shapes = tmpShape;
               } else {
                  ffpmsg( "Failed to allocate memory for Region data");
                  *status = MEMORY_ALLOCATION;
                  goto error;
               }

            }
            newShape        = &aRgn->Shapes[aRgn->nShapes++];
            newShape->sign  = 1;
            newShape->shape = point_rgn;

            /*  Check for format code at beginning of the line */

            if( !strncasecmp( namePtr, "image;", 6 ) ) {
				namePtr += 6;
				cFmt = pixel_fmt;
            } else if( !strncasecmp( namePtr, "fk4;", 4 ) ) {
				namePtr += 4;
				cFmt = degree_fmt;
            } else if( !strncasecmp( namePtr, "fk5;", 4 ) ) {
				namePtr += 4;
				cFmt = degree_fmt;
            } else if( !strncasecmp( namePtr, "icrs;", 5 ) ) {
				namePtr += 5;
				cFmt = degree_fmt;
            } else if( !strncasecmp( namePtr, "fk5", 3 ) ) {
				cFmt = degree_fmt;
                                continue;  /* supports POW region file format */
            } else if( !strncasecmp( namePtr, "fk4", 3 ) ) {
				cFmt = degree_fmt;
                                continue;  /* supports POW region file format */
            } else if( !strncasecmp( namePtr, "icrs", 4 ) ) {
				cFmt = degree_fmt;
                                continue;  /* supports POW region file format */
            } else if( !strncasecmp( namePtr, "image", 5 ) ) {
				cFmt = pixel_fmt;
                                continue;  /* supports POW region file format */
            } else if( !strncasecmp( namePtr, "galactic;", 9 ) ) {
               ffpmsg( "Galactic region coordinates not supported" );
               ffpmsg( namePtr );
			   *status = PARSE_SYNTAX_ERR;
               goto error;
			} else if( !strncasecmp( namePtr, "ecliptic;", 9 ) ) {
               ffpmsg( "ecliptic region coordinates not supported" );
               ffpmsg( namePtr );
			   *status = PARSE_SYNTAX_ERR;
               goto error;
            }

            while( *namePtr==' ' ) namePtr++;
            
			/*  Check for the shape's sign  */

            if( *namePtr=='+' ) {
               namePtr++;
            } else if( *namePtr=='-' ) {
               namePtr++;
               newShape->sign = 0;
            }

            /* Skip white space in region name */

            while( *namePtr==' ' ) namePtr++;
            if( *namePtr=='\0' ) {
               ffpmsg( "Error reading Region file" );
               *status = PARSE_SYNTAX_ERR;
               goto error;
            }
            lineLen = strlen( namePtr ) - 1;
            while( namePtr[lineLen]==' ' ) namePtr[lineLen--] = '\0';

            /*  Now identify the region  */

            if(        !strcasecmp( namePtr, "circle"  ) ) {
               newShape->shape = circle_rgn;
               if( nParams != 3 )
                  *status = PARSE_SYNTAX_ERR;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "annulus" ) ) {
               newShape->shape = annulus_rgn;
               if( nParams != 4 )
                  *status = PARSE_SYNTAX_ERR;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "ellipse" ) ) {
               newShape->shape = ellipse_rgn;
               if( nParams < 4 || nParams > 5 )
                  *status = PARSE_SYNTAX_ERR;
               newShape->param.gen.p[4] = 0.0;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "elliptannulus" ) ) {
               newShape->shape = elliptannulus_rgn;
               if( !( nParams==8 || nParams==6 ) )
                  *status = PARSE_SYNTAX_ERR;
               newShape->param.gen.p[6] = 0.0;
               newShape->param.gen.p[7] = 0.0;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "box"    ) 
                    || !strcasecmp( namePtr, "rotbox" ) ) {
               newShape->shape = box_rgn;
               if( nParams < 4 || nParams > 5 )
                  *status = PARSE_SYNTAX_ERR;
               newShape->param.gen.p[4] = 0.0;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "rectangle"    )
                    || !strcasecmp( namePtr, "rotrectangle" ) ) {
               newShape->shape = rectangle_rgn;
               if( nParams < 4 || nParams > 5 )
                  *status = PARSE_SYNTAX_ERR;
               newShape->param.gen.p[4] = 0.0;
               nCoords = 4;
            } else if( !strcasecmp( namePtr, "diamond"    )
                    || !strcasecmp( namePtr, "rotdiamond" )
                    || !strcasecmp( namePtr, "rhombus"    )
                    || !strcasecmp( namePtr, "rotrhombus" ) ) {
               newShape->shape = diamond_rgn;
               if( nParams < 4 || nParams > 5 )
                  *status = PARSE_SYNTAX_ERR;
               newShape->param.gen.p[4] = 0.0;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "sector"  )
                    || !strcasecmp( namePtr, "pie"     ) ) {
               newShape->shape = sector_rgn;
               if( nParams != 4 )
                  *status = PARSE_SYNTAX_ERR;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "point"   ) ) {
               newShape->shape = point_rgn;
               if( nParams != 2 )
                  *status = PARSE_SYNTAX_ERR;
               nCoords = 2;
            } else if( !strcasecmp( namePtr, "line"    ) ) {
               newShape->shape = line_rgn;
               if( nParams != 4 )
                  *status = PARSE_SYNTAX_ERR;
               nCoords = 4;
            } else if( !strcasecmp( namePtr, "polygon" ) ) {
               newShape->shape = poly_rgn;
               if( nParams < 6 || (nParams&1) )
                  *status = PARSE_SYNTAX_ERR;
               nCoords = nParams;
            } else {
               ffpmsg( "Unrecognized region found in region file:" );
               ffpmsg( namePtr );
               *status = PARSE_SYNTAX_ERR;
               goto error;
            }
            if( *status ) {
               ffpmsg( "Wrong number of parameters found for region" );
               ffpmsg( namePtr );
               goto error;
            }

            /*  Parse Parameter string... convert to pixels if necessary  */

            if( newShape->shape==poly_rgn ) {
               newShape->param.poly.Pts = (double *)malloc( nParams
                                                            * sizeof(double) );
               if( !newShape->param.poly.Pts ) {
                  ffpmsg(
                      "Could not allocate memory to hold polygon parameters" );
                  *status = MEMORY_ALLOCATION;
                  goto error;
               }
               newShape->param.poly.nPts = nParams;
               coords = newShape->param.poly.Pts;
            } else
               coords = newShape->param.gen.p;

            /*  Parse the initial "WCS?" coordinates  */
            for( i=0; i<nCoords; i+=2 ) {
               pX = paramPtr;
               while( *paramPtr!=',' ) paramPtr++;
               *(paramPtr++) = '\0';

               pY = paramPtr;
               while( *paramPtr!=',' && *paramPtr != '\0' ) paramPtr++;
               *(paramPtr++) = '\0';

               if( strchr(pX, ':' ) ) {
                  /*  Read in special format & convert to decimal degrees  */
                  cFmt = hhmmss_fmt;
                  mm = 0;
                  ss = 0.;
                  hh = strtol(pX, &endp, 10);
                  if (endp && *endp==':') {
                      pX = endp + 1;
                      mm = strtol(pX, &endp, 10);
                      if (endp && *endp==':') {
                          pX = endp + 1;
                          ss = atof( pX );
                      }
                  }
                  X = 15. * (hh + mm/60. + ss/3600.); /* convert to degrees */

                  mm = 0;
                  ss = 0.;
                  negdec = 0;

                  while( *pY==' ' ) pY++;
                  if (*pY=='-') {
                      negdec = 1;
                      pY++;
                  }
                  dd = strtol(pY, &endp, 10);
                  if (endp && *endp==':') {
                      pY = endp + 1;
                      mm = strtol(pY, &endp, 10);
                      if (endp && *endp==':') {
                          pY = endp + 1;
                          ss = atof( pY );
                      }
                  }
                  if (negdec)
                     Y = -dd - mm/60. - ss/3600.; /* convert to degrees */
                  else
                     Y = dd + mm/60. + ss/3600.;

               } else {
                  X = atof( pX );
                  Y = atof( pY );
               }
               if (i==0) {   /* save 1st coord. in case needed later */
                   xsave = X;
                   ysave = Y;
               }

               if( cFmt!=pixel_fmt ) {
                  /*  Convert to pixels  */
                  if( wcs==NULL || ! wcs->exists ) {
                     ffpmsg("WCS information needed to convert region coordinates.");
                     *status = NO_WCS_KEY;
                     goto error;
                  }
                  
                  if( ffxypx(  X,  Y, wcs->xrefval, wcs->yrefval,
                                      wcs->xrefpix, wcs->yrefpix,
                                      wcs->xinc,    wcs->yinc,
                                      wcs->rot,     wcs->type,
                              &x, &y, status ) ) {
                     ffpmsg("Error converting region to pixel coordinates.");
                     goto error;
                  }
                  X = x; Y = y;
               }
               coords[i]   = X;
               coords[i+1] = Y;
            }

            /*  Read in remaining parameters...  */

            for( ; i<nParams; i++ ) {
               pX = paramPtr;
               while( *paramPtr!=',' && *paramPtr != '\0' ) paramPtr++;
               *(paramPtr++) = '\0';
               coords[i] = strtod( pX, &endp );

               if (endp && *endp=='"') {
                   /* parameter given in arcsec so convert to pixels. */
                   /* Increment first Y coordinate by this amount then calc */
                   /* the distance in pixels from the original coordinate. */
                   /* NOTE: This assumes the pixels are square!! */
                   if (ysave < 0.)
                       Y = ysave + coords[i]/3600.;  /* don't exceed -90 */
                   else
                       Y = ysave - coords[i]/3600.;  /* don't exceed +90 */

                   X = xsave;
                   if( ffxypx(  X,  Y, wcs->xrefval, wcs->yrefval,
                                      wcs->xrefpix, wcs->yrefpix,
                                      wcs->xinc,    wcs->yinc,
                                      wcs->rot,     wcs->type,
                              &x, &y, status ) ) {
                     ffpmsg("Error converting region to pixel coordinates.");
                     goto error;
                   }

                   coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );

               } else if (endp && *endp=='\'') {
                   /* parameter given in arcmin so convert to pixels. */
                   /* Increment first Y coordinate by this amount, then calc */
                   /* the distance in pixels from the original coordinate. */
                   /* NOTE: This assumes the pixels are square!! */
                   if (ysave < 0.)
                       Y = ysave + coords[i]/60.;  /* don't exceed -90 */
                   else
                       Y = ysave - coords[i]/60.;  /* don't exceed +90 */

                   X = xsave;
                   if( ffxypx(  X,  Y, wcs->xrefval, wcs->yrefval,
                                      wcs->xrefpix, wcs->yrefpix,
                                      wcs->xinc,    wcs->yinc,
                                      wcs->rot,     wcs->type,
                              &x, &y, status ) ) {
                     ffpmsg("Error converting region to pixel coordinates.");
                     goto error;
                   }

                   coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );

               } else if (endp && *endp=='d') {
                   /* parameter given in degrees so convert to pixels. */
                   /* Increment first Y coordinate by this amount, then calc */
                   /* the distance in pixels from the original coordinate. */
                   /* NOTE: This assumes the pixels are square!! */
                   if (ysave < 0.)
                       Y = ysave + coords[i];  /* don't exceed -90 */
                   else
                       Y = ysave - coords[i];  /* don't exceed +90 */

                   X = xsave;
                   if( ffxypx(  X,  Y, wcs->xrefval, wcs->yrefval,
                                      wcs->xrefpix, wcs->yrefpix,
                                      wcs->xinc,    wcs->yinc,
                                      wcs->rot,     wcs->type,
                              &x, &y, status ) ) {
                     ffpmsg("Error converting region to pixel coordinates.");
                     goto error;
                   }

                   coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );
               }
            }

            /* Perform some useful calculations now to speed up filter later */

            switch( newShape->shape ) {
            case circle_rgn:
               newShape->param.gen.a = coords[2] * coords[2];
               break;
            case annulus_rgn:
               newShape->param.gen.a = coords[2] * coords[2];
               newShape->param.gen.b = coords[3] * coords[3];
               break;
            case sector_rgn:
               while( coords[2]> 180.0 ) coords[2] -= 360.0;
               while( coords[2]<=-180.0 ) coords[2] += 360.0;
               while( coords[3]> 180.0 ) coords[3] -= 360.0;
               while( coords[3]<=-180.0 ) coords[3] += 360.0;
               break;
            case ellipse_rgn:
               newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
               newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
               break;
            case elliptannulus_rgn:
               newShape->param.gen.a    = sin( myPI * (coords[6] / 180.0) );
               newShape->param.gen.b    = cos( myPI * (coords[6] / 180.0) );
               newShape->param.gen.sinT = sin( myPI * (coords[7] / 180.0) );
               newShape->param.gen.cosT = cos( myPI * (coords[7] / 180.0) );
               break;
            case box_rgn:
               newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
               newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
               break;
            case rectangle_rgn:
               newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
               newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
               X = 0.5 * ( coords[2]-coords[0] );
               Y = 0.5 * ( coords[3]-coords[1] );
               newShape->param.gen.a = fabs( X * newShape->param.gen.cosT
                                             + Y * newShape->param.gen.sinT );
               newShape->param.gen.b = fabs( Y * newShape->param.gen.cosT
                                             - X * newShape->param.gen.sinT );
               newShape->param.gen.p[5] = 0.5 * ( coords[2]+coords[0] );
               newShape->param.gen.p[6] = 0.5 * ( coords[3]+coords[1] );
               break;
            case diamond_rgn:
               newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
               newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
               break;
            case line_rgn:
               X = coords[2] - coords[0];
               Y = coords[3] - coords[1];
               R = sqrt( X*X + Y*Y );
               newShape->param.gen.sinT = ( R ? Y/R : 0.0 );
               newShape->param.gen.cosT = ( R ? X/R : 1.0 );
               newShape->param.gen.a    = R + 0.5;
               break;
            case point_rgn:
               break;
            case poly_rgn:
               /*  Find bounding box  */
               newShape->param.poly.xmin = coords[0];
               newShape->param.poly.xmax = coords[0];
               newShape->param.poly.ymin = coords[1];
               newShape->param.poly.ymax = coords[1];
               for( i=2; i<nParams; ) {
                  if( newShape->param.poly.xmin > coords[i] ) /* Min X */
                      newShape->param.poly.xmin = coords[i];
                  if( newShape->param.poly.xmax < coords[i] ) /* Max X */
                      newShape->param.poly.xmax = coords[i];
                  i++;
                  if( newShape->param.poly.ymin > coords[i] ) /* Min Y */
                      newShape->param.poly.ymin = coords[i];
                  if( newShape->param.poly.ymax < coords[i] ) /* Max Y */
                      newShape->param.poly.ymax = coords[i];
                  i++;
               }
               break;
            }

         }  /* End of while( *currLoc ) */
/*
  if (coords)printf("%.8f %.8f %.8f %.8f %.8f\n",
   coords[0],coords[1],coords[2],coords[3],coords[4]); 
*/
      }  /* End of if...else parse line */
   }   /* End of while( fgets(rgnFile) ) */


error:

   if( *status )
      fits_free_region( aRgn );
   else
      *Rgn = aRgn;

   fclose( rgnFile );
   free( currLine );

   return( *status );
}

/*---------------------------------------------------------------------------*/
int fftrgn( double    X,
            double    Y,
            SAORegion *Rgn )
/*  Test if the given point is within the region described by Rgn.  X and    */
/*  Y are in pixel coordinates.                                              */
/*---------------------------------------------------------------------------*/
{
   double x, y, dx, dy, xprime, yprime, r;
   RgnShape *Shapes;
   int i;
   int result = 0;

   Shapes = Rgn->Shapes;

   /* if an excluded region is given first, then implicitly   */
   /* assume a previous shape that includes the entire image. */
   if (!Shapes->sign)
      result = 1;

   for( i=0; i<Rgn->nShapes; i++, Shapes++ ) {

    /* only need to test if  */
    /*   the point is not already included and this is an include region, */
    /* or the point is included and this is an excluded region */

    if ( (!result && Shapes->sign) || (result && !Shapes->sign) ) { 

      result = 1;

      switch( Shapes->shape ) {

      case box_rgn:
         /*  Shift origin to center of region  */
         xprime = X - Shapes->param.gen.p[0];
         yprime = Y - Shapes->param.gen.p[1];

         /*  Rotate point to region's orientation  */
         x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
         y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;

         dx = 0.5 * Shapes->param.gen.p[2];
         dy = 0.5 * Shapes->param.gen.p[3];
         if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
            result = 0;
         break;

      case rectangle_rgn:
         /*  Shift origin to center of region  */
         xprime = X - Shapes->param.gen.p[5];
         yprime = Y - Shapes->param.gen.p[6];

         /*  Rotate point to region's orientation  */
         x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
         y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;

         dx = Shapes->param.gen.a;
         dy = Shapes->param.gen.b;
         if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
            result = 0;
         break;

      case diamond_rgn:
         /*  Shift origin to center of region  */
         xprime = X - Shapes->param.gen.p[0];
         yprime = Y - Shapes->param.gen.p[1];

         /*  Rotate point to region's orientation  */
         x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
         y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;

         dx = 0.5 * Shapes->param.gen.p[2];
         dy = 0.5 * Shapes->param.gen.p[3];
         r  = fabs(x/dx) + fabs(y/dy);
         if( r > 1 )
            result = 0;
         break;

      case circle_rgn:
         /*  Shift origin to center of region  */
         x = X - Shapes->param.gen.p[0];
         y = Y - Shapes->param.gen.p[1];

         r  = x*x + y*y;
         if ( r > Shapes->param.gen.a )
            result = 0;
         break;

      case annulus_rgn:
         /*  Shift origin to center of region  */
         x = X - Shapes->param.gen.p[0];
         y = Y - Shapes->param.gen.p[1];

         r = x*x + y*y;
         if ( r < Shapes->param.gen.a || r > Shapes->param.gen.b )
            result = 0;
         break;

      case sector_rgn:
         /*  Shift origin to center of region  */
         x = X - Shapes->param.gen.p[0];
         y = Y - Shapes->param.gen.p[1];

         if( x || y ) {
            r = atan2( y, x ) * 180.0 / myPI;
            if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) {
               if( r < Shapes->param.gen.p[2] || r > Shapes->param.gen.p[3] )
                  result = 0;
            } else {
               if( r < Shapes->param.gen.p[2] && r > Shapes->param.gen.p[3] )
                  result = 0;
            }
         }
         break;

      case ellipse_rgn:
         /*  Shift origin to center of region  */
         xprime = X - Shapes->param.gen.p[0];
         yprime = Y - Shapes->param.gen.p[1];

         /*  Rotate point to region's orientation  */
         x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
         y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;

         x /= Shapes->param.gen.p[2];
         y /= Shapes->param.gen.p[3];
         r = x*x + y*y;
         if( r>1.0 )
            result = 0;
         break;

      case elliptannulus_rgn:
         /*  Shift origin to center of region  */
         xprime = X - Shapes->param.gen.p[0];
         yprime = Y - Shapes->param.gen.p[1];

         /*  Rotate point to outer ellipse's orientation  */
         x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
         y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;

         x /= Shapes->param.gen.p[4];
         y /= Shapes->param.gen.p[5];
         r = x*x + y*y;
         if( r>1.0 )
            result = 0;
         else {
            /*  Repeat test for inner ellipse  */
            x =  xprime * Shapes->param.gen.b + yprime * Shapes->param.gen.a;
            y = -xprime * Shapes->param.gen.a + yprime * Shapes->param.gen.b;

            x /= Shapes->param.gen.p[2];
            y /= Shapes->param.gen.p[3];
            r = x*x + y*y;
            if( r<1.0 )
               result = 0;
         }
         break;

      case line_rgn:
         /*  Shift origin to first point of line  */
         xprime = X - Shapes->param.gen.p[0];
         yprime = Y - Shapes->param.gen.p[1];

         /*  Rotate point to line's orientation  */
         x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
         y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;

         if( (y < -0.5) || (y >= 0.5) || (x < -0.5)
             || (x >= Shapes->param.gen.a) )
            result = 0;
         break;

      case point_rgn:
         /*  Shift origin to center of region  */
         x = X - Shapes->param.gen.p[0];
         y = Y - Shapes->param.gen.p[1];

         if ( (x<-0.5) || (x>=0.5) || (y<-0.5) || (y>=0.5) )
            result = 0;
         break;

      case poly_rgn:
         if( X<Shapes->param.poly.xmin || X>Shapes->param.poly.xmax
             || Y<Shapes->param.poly.ymin || Y>Shapes->param.poly.ymax )
            result = 0;
         else
            result = Pt_in_Poly( X, Y, Shapes->param.poly.nPts,
                                       Shapes->param.poly.Pts );
         break;
      }

      if( !Shapes->sign ) result = !result;

     } 
   }

   return( result );
}

/*---------------------------------------------------------------------------*/
void fffrgn( SAORegion *Rgn )
/*   Free up memory allocated to hold the region data.                       */
/*---------------------------------------------------------------------------*/
{
   int i;

   for( i=0; i<Rgn->nShapes; i++ )
      if( Rgn->Shapes[i].shape == poly_rgn )
         free( Rgn->Shapes[i].param.poly.Pts );
   if( Rgn->Shapes )
      free( Rgn->Shapes );
   free( Rgn );
}

/*---------------------------------------------------------------------------*/
static int Pt_in_Poly( double x,
                       double y,
                       int nPts,
                       double *Pts )
/*  Internal routine for testing whether the coordinate x,y is within the    */
/*  polygon region traced out by the array Pts.                              */
/*---------------------------------------------------------------------------*/
{
   int i, j, flag=0;
   double prevX, prevY;
   double nextX, nextY;
   double dx, dy, Dy;

   nextX = Pts[nPts-2];
   nextY = Pts[nPts-1];

   for( i=0; i<nPts; i+=2 ) {
      prevX = nextX;
      prevY = nextY;

      nextX = Pts[i];
      nextY = Pts[i+1];

      if( (y>prevY && y>=nextY) || (y<prevY && y<=nextY)
          || (x>prevX && x>=nextX) )
         continue;
      
      /* Check to see if x,y lies right on the segment */

      if( x>=prevX || x>nextX ) {
         dy = y - prevY;
         Dy = nextY - prevY;

         if( fabs(Dy)<1e-10 ) {
            if( fabs(dy)<1e-10 )
               return( 1 );
            else
               continue;
         }

         dx = prevX + ( (nextX-prevX)/(Dy) ) * dy - x;
         if( dx < -1e-10 )
            continue;
         if( dx <  1e-10 )
            return( 1 );
      }

      /* There is an intersection! Make sure it isn't a V point.  */

      if( y != prevY ) {
         flag = 1 - flag;
      } else {
         j = i+1;  /* Point to Y component */
         do {
            if( j>1 )
               j -= 2;
            else
               j = nPts-1;
         } while( y == Pts[j] );

         if( (nextY-y)*(y-Pts[j]) > 0 )
            flag = 1-flag;
      }

   }
   return( flag );
}