Initial commit

1 files changed, 919 insertions, 0 deletions

diff --git a/pkg/tbtables/cfitsio/region.c b/pkg/tbtables/cfitsio/region.c
new file mode 100644
index 00000000..31df9b7a
--- /dev/null
+++ b/pkg/tbtables/cfitsio/region.c
@@ -0,0 +1,919 @@
+#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;
+
+
+            /*  Check for format code at beginning of the line */
+
+            if( !strncasecmp( namePtr, "image;", 6 ) ) {
+				namePtr += 6;
+				cFmt = pixel_fmt;
+            } else if( !strncasecmp( namePtr, "physical;", 9 ) ) {
+                                namePtr += 9;
+                                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;
+
+            /* the following 4 cases support region files created by POW which
+               may have lines containing  only a format code, not followed
+               by a ';' (and with no region specifier on the line).  We use
+               the 'continue' statement to jump to the end of the loop and
+               then continue reading the next line of the region file. */
+
+            } 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;
+            }
+
+            /**************************************************/
+            /*  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;
+
+            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 );
+}
+