"Fossies" - the Fresh Open Source Software Archive

Member "cfitsio-4.0.0/region.c" (13 Jul 2021, 57375 Bytes) of package /linux/misc/cfitsio-4.0.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 #include <stdio.h>
    2 #include <stdlib.h>
    3 #include <string.h>
    4 #include <math.h>
    5 #include <ctype.h>
    6 #include "fitsio2.h"
    7 #include "region.h"
    8 static int Pt_in_Poly( double x, double y, int nPts, double *Pts );
    9 
   10 /*---------------------------------------------------------------------------*/
   11 int fits_read_rgnfile( const char *filename,
   12             WCSdata    *wcs,
   13             SAORegion  **Rgn,
   14             int        *status )
   15 /*  Read regions from either a FITS or ASCII region file and return the information     */
   16 /*  in the "SAORegion" structure.  If it is nonNULL, use wcs to convert the  */
   17 /*  region coordinates to pixels.  Return an error if region is in degrees   */
   18 /*  but no WCS data is provided.                                             */
   19 /*---------------------------------------------------------------------------*/
   20 {
   21   fitsfile *fptr;
   22   int tstatus = 0;
   23 
   24   if( *status ) return( *status );
   25 
   26   /* try to open as a FITS file - if that doesn't work treat as an ASCII file */
   27 
   28   fits_write_errmark();
   29   if ( ffopen(&fptr, filename, READONLY, &tstatus) ) {
   30     fits_clear_errmark();
   31     fits_read_ascii_region(filename, wcs, Rgn, status);
   32   } else {
   33     fits_read_fits_region(fptr, wcs, Rgn, status);
   34   }
   35 
   36   return(*status);
   37 
   38 }
   39 /*---------------------------------------------------------------------------*/
   40 int fits_read_ascii_region( const char *filename,
   41                 WCSdata    *wcs,
   42                 SAORegion  **Rgn,
   43                 int        *status )
   44 /*  Read regions from a SAO-style region file and return the information     */
   45 /*  in the "SAORegion" structure.  If it is nonNULL, use wcs to convert the  */
   46 /*  region coordinates to pixels.  Return an error if region is in degrees   */
   47 /*  but no WCS data is provided.                                             */
   48 /*---------------------------------------------------------------------------*/
   49 {
   50    char     *currLine;
   51    char     *namePtr, *paramPtr, *currLoc;
   52    char     *pX, *pY, *endp;
   53    long     allocLen, lineLen, hh, mm, dd;
   54    double   *coords, X, Y, x, y, ss, div, xsave= 0., ysave= 0.;
   55    int      nParams, nCoords, negdec;
   56    int      i, done;
   57    FILE     *rgnFile;
   58    coordFmt cFmt;
   59    SAORegion *aRgn;
   60    RgnShape *newShape, *tmpShape;
   61 
   62    if( *status ) return( *status );
   63 
   64    aRgn = (SAORegion *)malloc( sizeof(SAORegion) );
   65    if( ! aRgn ) {
   66       ffpmsg("Couldn't allocate memory to hold Region file contents.");
   67       return(*status = MEMORY_ALLOCATION );
   68    }
   69    aRgn->nShapes    =    0;
   70    aRgn->Shapes     = NULL;
   71    if( wcs && wcs->exists )
   72       aRgn->wcs = *wcs;
   73    else
   74       aRgn->wcs.exists = 0;
   75 
   76    cFmt = pixel_fmt; /* set default format */
   77 
   78    /*  Allocate Line Buffer  */
   79 
   80    allocLen = 512;
   81    currLine = (char *)malloc( allocLen * sizeof(char) );
   82    if( !currLine ) {
   83       free( aRgn );
   84       ffpmsg("Couldn't allocate memory to hold Region file contents.");
   85       return(*status = MEMORY_ALLOCATION );
   86    }
   87 
   88    /*  Open Region File  */
   89 
   90    if( (rgnFile = fopen( filename, "r" ))==NULL ) {
   91       snprintf(currLine,allocLen,"Could not open Region file %s.",filename);
   92       ffpmsg( currLine );
   93       free( currLine );
   94       free( aRgn );
   95       return( *status = FILE_NOT_OPENED );
   96    }
   97    
   98    /*  Read in file, line by line  */
   99    /*  First, set error status in case file is empty */ 
  100    *status = FILE_NOT_OPENED;
  101 
  102    while( fgets(currLine,allocLen,rgnFile) != NULL ) {
  103 
  104       /* reset status if we got here */
  105       *status = 0;
  106 
  107       /*  Make sure we have a full line of text  */
  108 
  109       lineLen = strlen(currLine);
  110       while( lineLen==allocLen-1 && currLine[lineLen-1]!='\n' ) {
  111          currLoc = (char *)realloc( currLine, 2 * allocLen * sizeof(char) );
  112          if( !currLoc ) {
  113             ffpmsg("Couldn't allocate memory to hold Region file contents.");
  114             *status = MEMORY_ALLOCATION;
  115             goto error;
  116          } else {
  117             currLine = currLoc;
  118          }
  119          fgets( currLine+lineLen, allocLen+1, rgnFile );
  120          allocLen += allocLen;
  121          lineLen  += strlen(currLine+lineLen);
  122       }
  123 
  124       currLoc = currLine;
  125       if( *currLoc == '#' ) {
  126 
  127          /*  Look to see if it is followed by a format statement...  */
  128          /*  if not skip line                                        */
  129 
  130          currLoc++;
  131          while( isspace(*currLoc) ) currLoc++;
  132          if( !fits_strncasecmp( currLoc, "format:", 7 ) ) {
  133             if( aRgn->nShapes ) {
  134                ffpmsg("Format code encountered after reading 1 or more shapes.");
  135                *status = PARSE_SYNTAX_ERR;
  136                goto error;
  137             }
  138             currLoc += 7;
  139             while( isspace(*currLoc) ) currLoc++;
  140             if( !fits_strncasecmp( currLoc, "pixel", 5 ) ) {
  141                cFmt = pixel_fmt;
  142             } else if( !fits_strncasecmp( currLoc, "degree", 6 ) ) {
  143                cFmt = degree_fmt;
  144             } else if( !fits_strncasecmp( currLoc, "hhmmss", 6 ) ) {
  145                cFmt = hhmmss_fmt;
  146             } else if( !fits_strncasecmp( currLoc, "hms", 3 ) ) {
  147                cFmt = hhmmss_fmt;
  148             } else {
  149                ffpmsg("Unknown format code encountered in region file.");
  150                *status = PARSE_SYNTAX_ERR;
  151                goto error;
  152             }
  153          }
  154 
  155       } else if( !fits_strncasecmp( currLoc, "glob", 4 ) ) {
  156           /* skip lines that begin with the word 'global' */
  157 
  158       } else {
  159 
  160          while( *currLoc != '\0' ) {
  161 
  162             namePtr  = currLoc;
  163             paramPtr = NULL;
  164             nParams  = 1;
  165 
  166             /*  Search for closing parenthesis  */
  167 
  168             done = 0;
  169             while( !done && !*status && *currLoc ) {
  170                switch (*currLoc) {
  171                case '(':
  172                   *currLoc = '\0';
  173                   currLoc++;
  174                   if( paramPtr )   /* Can't have two '(' in a region! */
  175                      *status = 1;
  176                   else
  177                      paramPtr = currLoc;
  178                   break;
  179                case ')':
  180                   *currLoc = '\0';
  181                   currLoc++;
  182                   if( !paramPtr )  /* Can't have a ')' without a '(' first */
  183                      *status = 1;
  184                   else
  185                      done = 1;
  186                   break;
  187                case '#':
  188                case '\n':
  189                   *currLoc = '\0';
  190                   if( !paramPtr )  /* Allow for a blank line */
  191                      done = 1;
  192                   break;
  193                case ':':  
  194                   currLoc++;
  195                   if ( paramPtr ) cFmt = hhmmss_fmt; /* set format if parameter has : */
  196                   break;
  197                case 'd':
  198                   currLoc++;
  199                   if ( paramPtr ) cFmt = degree_fmt; /* set format if parameter has d */  
  200                   break;
  201                case ',':
  202                   nParams++;  /* Fall through to default */
  203                default:
  204                   currLoc++;
  205                   break;
  206                }
  207             }
  208             if( *status || !done ) {
  209                ffpmsg( "Error reading Region file" );
  210                *status = PARSE_SYNTAX_ERR;
  211                goto error;
  212             }
  213 
  214             /*  Skip white space in region name  */
  215 
  216             while( isspace(*namePtr) ) namePtr++;
  217 
  218             /*  Was this a blank line? Or the end of the current one  */
  219 
  220             if( ! *namePtr && ! paramPtr ) continue;
  221 
  222             /*  Check for format code at beginning of the line */
  223 
  224             if( !fits_strncasecmp( namePtr, "image;", 6 ) ) {
  225                 namePtr += 6;
  226                 cFmt = pixel_fmt;
  227             } else if( !fits_strncasecmp( namePtr, "physical;", 9 ) ) {
  228                                 namePtr += 9;
  229                                 cFmt = pixel_fmt;
  230             } else if( !fits_strncasecmp( namePtr, "linear;", 7 ) ) {
  231                                 namePtr += 7;
  232                                 cFmt = pixel_fmt;
  233             } else if( !fits_strncasecmp( namePtr, "fk4;", 4 ) ) {
  234                 namePtr += 4;
  235                 cFmt = degree_fmt;
  236             } else if( !fits_strncasecmp( namePtr, "fk5;", 4 ) ) {
  237                 namePtr += 4;
  238                 cFmt = degree_fmt;
  239             } else if( !fits_strncasecmp( namePtr, "icrs;", 5 ) ) {
  240                 namePtr += 5;
  241                 cFmt = degree_fmt;
  242 
  243             /* the following 5 cases support region files created by POW 
  244            (or ds9 Version 4.x) which
  245                may have lines containing  only a format code, not followed
  246                by a ';' (and with no region specifier on the line).  We use
  247                the 'continue' statement to jump to the end of the loop and
  248                then continue reading the next line of the region file. */
  249 
  250             } else if( !fits_strncasecmp( namePtr, "fk5", 3 ) ) {
  251                 cFmt = degree_fmt;
  252                                 continue;  /* supports POW region file format */
  253             } else if( !fits_strncasecmp( namePtr, "fk4", 3 ) ) {
  254                 cFmt = degree_fmt;
  255                                 continue;  /* supports POW region file format */
  256             } else if( !fits_strncasecmp( namePtr, "icrs", 4 ) ) {
  257                 cFmt = degree_fmt;
  258                                 continue;  /* supports POW region file format */
  259             } else if( !fits_strncasecmp( namePtr, "image", 5 ) ) {
  260                 cFmt = pixel_fmt;
  261                                 continue;  /* supports POW region file format */
  262             } else if( !fits_strncasecmp( namePtr, "physical", 8 ) ) {
  263                 cFmt = pixel_fmt;
  264                                 continue;  /* supports POW region file format */
  265 
  266 
  267             } else if( !fits_strncasecmp( namePtr, "galactic;", 9 ) ) {
  268                ffpmsg( "Galactic region coordinates not supported" );
  269                ffpmsg( namePtr );
  270                *status = PARSE_SYNTAX_ERR;
  271                goto error;
  272             } else if( !fits_strncasecmp( namePtr, "ecliptic;", 9 ) ) {
  273                ffpmsg( "ecliptic region coordinates not supported" );
  274                ffpmsg( namePtr );
  275                *status = PARSE_SYNTAX_ERR;
  276                goto error;
  277             }
  278 
  279             /**************************************************/
  280             /*  We've apparently found a region... Set it up  */
  281             /**************************************************/
  282 
  283             if( !(aRgn->nShapes % 10) ) {
  284                if( aRgn->Shapes )
  285                   tmpShape = (RgnShape *)realloc( aRgn->Shapes,
  286                                                   (10+aRgn->nShapes)
  287                                                   * sizeof(RgnShape) );
  288                else
  289                   tmpShape = (RgnShape *) malloc( 10 * sizeof(RgnShape) );
  290                if( tmpShape ) {
  291                   aRgn->Shapes = tmpShape;
  292                } else {
  293                   ffpmsg( "Failed to allocate memory for Region data");
  294                   *status = MEMORY_ALLOCATION;
  295                   goto error;
  296                }
  297 
  298             }
  299             newShape        = &aRgn->Shapes[aRgn->nShapes++];
  300             newShape->sign  = 1;
  301             newShape->shape = point_rgn;
  302         for (i=0; i<8; i++) newShape->param.gen.p[i] = 0.0;
  303         newShape->param.gen.a = 0.0;
  304         newShape->param.gen.b = 0.0;
  305         newShape->param.gen.sinT = 0.0;
  306         newShape->param.gen.cosT = 0.0;
  307 
  308             while( isspace(*namePtr) ) namePtr++;
  309             
  310             /*  Check for the shape's sign  */
  311 
  312             if( *namePtr=='+' ) {
  313                namePtr++;
  314             } else if( *namePtr=='-' ) {
  315                namePtr++;
  316                newShape->sign = 0;
  317             }
  318 
  319             /* Skip white space in region name */
  320 
  321             while( isspace(*namePtr) ) namePtr++;
  322             if( *namePtr=='\0' ) {
  323                ffpmsg( "Error reading Region file" );
  324                *status = PARSE_SYNTAX_ERR;
  325                goto error;
  326             }
  327             lineLen = strlen( namePtr ) - 1;
  328             while( isspace(namePtr[lineLen]) ) namePtr[lineLen--] = '\0';
  329 
  330             /*  Now identify the region  */
  331 
  332             if(        !fits_strcasecmp( namePtr, "circle"  ) ) {
  333                newShape->shape = circle_rgn;
  334                if( nParams != 3 )
  335                   *status = PARSE_SYNTAX_ERR;
  336                nCoords = 2;
  337             } else if( !fits_strcasecmp( namePtr, "annulus" ) ) {
  338                newShape->shape = annulus_rgn;
  339                if( nParams != 4 )
  340                   *status = PARSE_SYNTAX_ERR;
  341                nCoords = 2;
  342             } else if( !fits_strcasecmp( namePtr, "ellipse" ) ) {
  343                if( nParams < 4 || nParams > 8 ) {
  344                   *status = PARSE_SYNTAX_ERR;
  345            } else if ( nParams < 6 ) {
  346          newShape->shape = ellipse_rgn;
  347          newShape->param.gen.p[4] = 0.0;
  348            } else {
  349          newShape->shape = elliptannulus_rgn;
  350          newShape->param.gen.p[6] = 0.0;
  351          newShape->param.gen.p[7] = 0.0;
  352            }
  353                nCoords = 2;
  354             } else if( !fits_strcasecmp( namePtr, "elliptannulus" ) ) {
  355                newShape->shape = elliptannulus_rgn;
  356                if( !( nParams==8 || nParams==6 ) )
  357                   *status = PARSE_SYNTAX_ERR;
  358                newShape->param.gen.p[6] = 0.0;
  359                newShape->param.gen.p[7] = 0.0;
  360                nCoords = 2;
  361             } else if( !fits_strcasecmp( namePtr, "box"    ) 
  362                     || !fits_strcasecmp( namePtr, "rotbox" ) ) {
  363            if( nParams < 4 || nParams > 8 ) {
  364          *status = PARSE_SYNTAX_ERR;
  365            } else if ( nParams < 6 ) {
  366          newShape->shape = box_rgn;
  367          newShape->param.gen.p[4] = 0.0;
  368            } else {
  369           newShape->shape = boxannulus_rgn;
  370           newShape->param.gen.p[6] = 0.0;
  371           newShape->param.gen.p[7] = 0.0;
  372            }
  373            nCoords = 2;
  374             } else if( !fits_strcasecmp( namePtr, "rectangle"    )
  375                     || !fits_strcasecmp( namePtr, "rotrectangle" ) ) {
  376                newShape->shape = rectangle_rgn;
  377                if( nParams < 4 || nParams > 5 )
  378                   *status = PARSE_SYNTAX_ERR;
  379                newShape->param.gen.p[4] = 0.0;
  380                nCoords = 4;
  381             } else if( !fits_strcasecmp( namePtr, "diamond"    )
  382                     || !fits_strcasecmp( namePtr, "rotdiamond" )
  383                     || !fits_strcasecmp( namePtr, "rhombus"    )
  384                     || !fits_strcasecmp( namePtr, "rotrhombus" ) ) {
  385                newShape->shape = diamond_rgn;
  386                if( nParams < 4 || nParams > 5 )
  387                   *status = PARSE_SYNTAX_ERR;
  388                newShape->param.gen.p[4] = 0.0;
  389                nCoords = 2;
  390             } else if( !fits_strcasecmp( namePtr, "sector"  )
  391                     || !fits_strcasecmp( namePtr, "pie"     ) ) {
  392                newShape->shape = sector_rgn;
  393                if( nParams != 4 )
  394                   *status = PARSE_SYNTAX_ERR;
  395                nCoords = 2;
  396             } else if( !fits_strcasecmp( namePtr, "point"   ) ) {
  397                newShape->shape = point_rgn;
  398                if( nParams != 2 )
  399                   *status = PARSE_SYNTAX_ERR;
  400                nCoords = 2;
  401             } else if( !fits_strcasecmp( namePtr, "line"    ) ) {
  402                newShape->shape = line_rgn;
  403                if( nParams != 4 )
  404                   *status = PARSE_SYNTAX_ERR;
  405                nCoords = 4;
  406             } else if( !fits_strcasecmp( namePtr, "polygon" ) ) {
  407                newShape->shape = poly_rgn;
  408                if( nParams < 6 || (nParams&1) )
  409                   *status = PARSE_SYNTAX_ERR;
  410                nCoords = nParams;
  411             } else if( !fits_strcasecmp( namePtr, "panda" ) ) {
  412                newShape->shape = panda_rgn;
  413                if( nParams != 8 )
  414                   *status = PARSE_SYNTAX_ERR;
  415                nCoords = 2;
  416             } else if( !fits_strcasecmp( namePtr, "epanda" ) ) {
  417                newShape->shape = epanda_rgn;
  418                if( nParams < 10 || nParams > 11 )
  419                   *status = PARSE_SYNTAX_ERR;
  420                newShape->param.gen.p[10] = 0.0;
  421                nCoords = 2;
  422             } else if( !fits_strcasecmp( namePtr, "bpanda" ) ) {
  423                newShape->shape = bpanda_rgn;
  424                if( nParams < 10 || nParams > 11 )
  425                   *status = PARSE_SYNTAX_ERR;
  426                newShape->param.gen.p[10] = 0.0;
  427                nCoords = 2;
  428             } else {
  429                ffpmsg( "Unrecognized region found in region file:" );
  430                ffpmsg( namePtr );
  431                *status = PARSE_SYNTAX_ERR;
  432                goto error;
  433             }
  434             if( *status ) {
  435                ffpmsg( "Wrong number of parameters found for region" );
  436                ffpmsg( namePtr );
  437                goto error;
  438             }
  439 
  440             /*  Parse Parameter string... convert to pixels if necessary  */
  441 
  442             if( newShape->shape==poly_rgn ) {
  443                newShape->param.poly.Pts = (double *)malloc( nParams
  444                                                             * sizeof(double) );
  445                if( !newShape->param.poly.Pts ) {
  446                   ffpmsg(
  447                       "Could not allocate memory to hold polygon parameters" );
  448                   *status = MEMORY_ALLOCATION;
  449                   goto error;
  450                }
  451                newShape->param.poly.nPts = nParams;
  452                coords = newShape->param.poly.Pts;
  453             } else
  454                coords = newShape->param.gen.p;
  455 
  456             /*  Parse the initial "WCS?" coordinates  */
  457             for( i=0; i<nCoords; i+=2 ) {
  458 
  459                pX = paramPtr;
  460                while( *paramPtr!=',' ) paramPtr++;
  461                *(paramPtr++) = '\0';
  462 
  463                pY = paramPtr;
  464                while( *paramPtr!=',' && *paramPtr != '\0' ) paramPtr++;
  465                *(paramPtr++) = '\0';
  466 
  467                if( strchr(pX, ':' ) ) {
  468                   /*  Read in special format & convert to decimal degrees  */
  469                   cFmt = hhmmss_fmt;
  470                   mm = 0;
  471                   ss = 0.;
  472                   hh = strtol(pX, &endp, 10);
  473                   if (endp && *endp==':') {
  474                       pX = endp + 1;
  475                       mm = strtol(pX, &endp, 10);
  476                       if (endp && *endp==':') {
  477                           pX = endp + 1;
  478                           ss = atof( pX );
  479                       }
  480                   }
  481                   X = 15. * (hh + mm/60. + ss/3600.); /* convert to degrees */
  482 
  483                   mm = 0;
  484                   ss = 0.;
  485                   negdec = 0;
  486 
  487                   while( isspace(*pY) ) pY++;
  488                   if (*pY=='-') {
  489                       negdec = 1;
  490                       pY++;
  491                   }
  492                   dd = strtol(pY, &endp, 10);
  493                   if (endp && *endp==':') {
  494                       pY = endp + 1;
  495                       mm = strtol(pY, &endp, 10);
  496                       if (endp && *endp==':') {
  497                           pY = endp + 1;
  498                           ss = atof( pY );
  499                       }
  500                   }
  501                   if (negdec)
  502                      Y = -dd - mm/60. - ss/3600.; /* convert to degrees */
  503                   else
  504                      Y = dd + mm/60. + ss/3600.;
  505 
  506                } else {
  507                   X = atof( pX );
  508                   Y = atof( pY );
  509                }
  510                if (i==0) {   /* save 1st coord. in case needed later */
  511                    xsave = X;
  512                    ysave = Y;
  513                }
  514 
  515                if( cFmt!=pixel_fmt ) {
  516                   /*  Convert to pixels  */
  517                   if( wcs==NULL || ! wcs->exists ) {
  518                      ffpmsg("WCS information needed to convert region coordinates.");
  519                      *status = NO_WCS_KEY;
  520                      goto error;
  521                   }
  522                   
  523                   if( ffxypx(  X,  Y, wcs->xrefval, wcs->yrefval,
  524                                       wcs->xrefpix, wcs->yrefpix,
  525                                       wcs->xinc,    wcs->yinc,
  526                                       wcs->rot,     wcs->type,
  527                               &x, &y, status ) ) {
  528                      ffpmsg("Error converting region to pixel coordinates.");
  529                      goto error;
  530                   }
  531                   X = x; Y = y;
  532                }
  533                coords[i]   = X;
  534                coords[i+1] = Y;
  535 
  536             }
  537 
  538             /*  Read in remaining parameters...  */
  539 
  540             for( ; i<nParams; i++ ) {
  541                pX = paramPtr;
  542                while( *paramPtr!=',' && *paramPtr != '\0' ) paramPtr++;
  543                *(paramPtr++) = '\0';
  544                coords[i] = strtod( pX, &endp );
  545 
  546            if (endp && (*endp=='"' || *endp=='\'' || *endp=='d') ) {
  547           div = 1.0;
  548           if ( *endp=='"' ) div = 3600.0;
  549           if ( *endp=='\'' ) div = 60.0;
  550           /* parameter given in arcsec so convert to pixels. */
  551           /* Increment first Y coordinate by this amount then calc */
  552           /* the distance in pixels from the original coordinate. */
  553           /* NOTE: This assumes the pixels are square!! */
  554           if (ysave < 0.)
  555              Y = ysave + coords[i]/div;  /* don't exceed -90 */
  556           else
  557              Y = ysave - coords[i]/div;  /* don't exceed +90 */
  558 
  559           X = xsave;
  560           if( ffxypx(  X,  Y, wcs->xrefval, wcs->yrefval,
  561                    wcs->xrefpix, wcs->yrefpix,
  562                    wcs->xinc,    wcs->yinc,
  563                    wcs->rot,     wcs->type,
  564                                &x, &y, status ) ) {
  565              ffpmsg("Error converting region to pixel coordinates.");
  566              goto error;
  567           }
  568          
  569           coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );
  570 
  571                }
  572             }
  573 
  574         /* special case for elliptannulus and boxannulus if only one angle
  575            was given */
  576 
  577         if ( (newShape->shape == elliptannulus_rgn || 
  578           newShape->shape == boxannulus_rgn ) && nParams == 7 ) {
  579           coords[7] = coords[6];
  580         }
  581 
  582             /* Also, correct the position angle for any WCS rotation:  */
  583             /*    If regions are specified in WCS coordintes, then the angles */
  584             /*    are relative to the WCS system, not the pixel X,Y system */
  585 
  586         if( cFmt!=pixel_fmt ) {     
  587           switch( newShape->shape ) {
  588           case sector_rgn:
  589           case panda_rgn:
  590         coords[2] += (wcs->rot);
  591         coords[3] += (wcs->rot);
  592         break;
  593           case box_rgn:
  594           case rectangle_rgn:
  595           case diamond_rgn:
  596           case ellipse_rgn:
  597         coords[4] += (wcs->rot);
  598         break;
  599           case boxannulus_rgn:
  600           case elliptannulus_rgn:
  601         coords[6] += (wcs->rot);
  602         coords[7] += (wcs->rot);
  603         break;
  604           case epanda_rgn:
  605           case bpanda_rgn:
  606         coords[2] += (wcs->rot);
  607         coords[3] += (wcs->rot);
  608         coords[10] += (wcs->rot);
  609               default:
  610                 break;
  611           }
  612         }
  613 
  614         /* do some precalculations to speed up tests */
  615 
  616         fits_setup_shape(newShape);
  617 
  618          }  /* End of while( *currLoc ) */
  619 /*
  620   if (coords)printf("%.8f %.8f %.8f %.8f %.8f\n",
  621    coords[0],coords[1],coords[2],coords[3],coords[4]); 
  622 */
  623       }  /* End of if...else parse line */
  624    }   /* End of while( fgets(rgnFile) ) */
  625 
  626    /* set up component numbers */
  627 
  628    fits_set_region_components( aRgn );
  629 
  630 error:
  631 
  632    if( *status ) {
  633       fits_free_region( aRgn );
  634    } else {
  635       *Rgn = aRgn;
  636    }
  637 
  638    fclose( rgnFile );
  639    free( currLine );
  640 
  641    return( *status );
  642 }
  643 
  644 /*---------------------------------------------------------------------------*/
  645 int fits_in_region( double    X,
  646             double    Y,
  647             SAORegion *Rgn )
  648 /*  Test if the given point is within the region described by Rgn.  X and    */
  649 /*  Y are in pixel coordinates.                                              */
  650 /*---------------------------------------------------------------------------*/
  651 {
  652    double x, y, dx, dy, xprime, yprime, r, th;
  653    RgnShape *Shapes;
  654    int i, cur_comp;
  655    int result, comp_result;
  656 
  657    Shapes = Rgn->Shapes;
  658 
  659    result = 0;
  660    comp_result = 0;
  661    cur_comp = Rgn->Shapes[0].comp;
  662 
  663    for( i=0; i<Rgn->nShapes; i++, Shapes++ ) {
  664 
  665      /* if this region has a different component number to the last one  */
  666      /* then replace the accumulated selection logical with the union of */
  667      /* the current logical and the total logical. Reinitialize the      */
  668      /* temporary logical.                                               */
  669 
  670      if ( i==0 || Shapes->comp != cur_comp ) {
  671        result = result || comp_result;
  672        cur_comp = Shapes->comp;
  673        /* if an excluded region is given first, then implicitly   */
  674        /* assume a previous shape that includes the entire image. */
  675        comp_result = !Shapes->sign;
  676      }
  677 
  678     /* only need to test if  */
  679     /*   the point is not already included and this is an include region, */
  680     /* or the point is included and this is an excluded region */
  681 
  682     if ( (!comp_result && Shapes->sign) || (comp_result && !Shapes->sign) ) { 
  683 
  684       comp_result = 1;
  685 
  686       switch( Shapes->shape ) {
  687 
  688       case box_rgn:
  689          /*  Shift origin to center of region  */
  690          xprime = X - Shapes->param.gen.p[0];
  691          yprime = Y - Shapes->param.gen.p[1];
  692 
  693          /*  Rotate point to region's orientation  */
  694          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  695          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  696 
  697          dx = 0.5 * Shapes->param.gen.p[2];
  698          dy = 0.5 * Shapes->param.gen.p[3];
  699          if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
  700             comp_result = 0;
  701          break;
  702 
  703       case boxannulus_rgn:
  704          /*  Shift origin to center of region  */
  705          xprime = X - Shapes->param.gen.p[0];
  706          yprime = Y - Shapes->param.gen.p[1];
  707 
  708          /*  Rotate point to region's orientation  */
  709          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  710          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  711 
  712          dx = 0.5 * Shapes->param.gen.p[4];
  713          dy = 0.5 * Shapes->param.gen.p[5];
  714          if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) ) {
  715        comp_result = 0;
  716      } else {
  717        /* Repeat test for inner box */
  718        x =  xprime * Shapes->param.gen.b + yprime * Shapes->param.gen.a;
  719        y = -xprime * Shapes->param.gen.a + yprime * Shapes->param.gen.b;
  720        
  721        dx = 0.5 * Shapes->param.gen.p[2];
  722        dy = 0.5 * Shapes->param.gen.p[3];
  723        if( (x >= -dx) && (x <= dx) && (y >= -dy) && (y <= dy) )
  724          comp_result = 0;
  725      }
  726          break;
  727 
  728       case rectangle_rgn:
  729          /*  Shift origin to center of region  */
  730          xprime = X - Shapes->param.gen.p[5];
  731          yprime = Y - Shapes->param.gen.p[6];
  732 
  733          /*  Rotate point to region's orientation  */
  734          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  735          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  736 
  737          dx = Shapes->param.gen.a;
  738          dy = Shapes->param.gen.b;
  739          if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
  740             comp_result = 0;
  741          break;
  742 
  743       case diamond_rgn:
  744          /*  Shift origin to center of region  */
  745          xprime = X - Shapes->param.gen.p[0];
  746          yprime = Y - Shapes->param.gen.p[1];
  747 
  748          /*  Rotate point to region's orientation  */
  749          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  750          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  751 
  752          dx = 0.5 * Shapes->param.gen.p[2];
  753          dy = 0.5 * Shapes->param.gen.p[3];
  754          r  = fabs(x/dx) + fabs(y/dy);
  755          if( r > 1 )
  756             comp_result = 0;
  757          break;
  758 
  759       case circle_rgn:
  760          /*  Shift origin to center of region  */
  761          x = X - Shapes->param.gen.p[0];
  762          y = Y - Shapes->param.gen.p[1];
  763 
  764          r  = x*x + y*y;
  765          if ( r > Shapes->param.gen.a )
  766             comp_result = 0;
  767          break;
  768 
  769       case annulus_rgn:
  770          /*  Shift origin to center of region  */
  771          x = X - Shapes->param.gen.p[0];
  772          y = Y - Shapes->param.gen.p[1];
  773 
  774          r = x*x + y*y;
  775          if ( r < Shapes->param.gen.a || r > Shapes->param.gen.b )
  776             comp_result = 0;
  777          break;
  778 
  779       case sector_rgn:
  780          /*  Shift origin to center of region  */
  781          x = X - Shapes->param.gen.p[0];
  782          y = Y - Shapes->param.gen.p[1];
  783 
  784          if( x || y ) {
  785             r = atan2( y, x ) * RadToDeg;
  786             if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) {
  787                if( r < Shapes->param.gen.p[2] || r > Shapes->param.gen.p[3] )
  788                   comp_result = 0;
  789             } else {
  790                if( r < Shapes->param.gen.p[2] && r > Shapes->param.gen.p[3] )
  791                   comp_result = 0;
  792             }
  793          }
  794          break;
  795 
  796       case ellipse_rgn:
  797          /*  Shift origin to center of region  */
  798          xprime = X - Shapes->param.gen.p[0];
  799          yprime = Y - Shapes->param.gen.p[1];
  800 
  801          /*  Rotate point to region's orientation  */
  802          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  803          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  804 
  805          x /= Shapes->param.gen.p[2];
  806          y /= Shapes->param.gen.p[3];
  807          r = x*x + y*y;
  808          if( r>1.0 )
  809             comp_result = 0;
  810          break;
  811 
  812       case elliptannulus_rgn:
  813          /*  Shift origin to center of region  */
  814          xprime = X - Shapes->param.gen.p[0];
  815          yprime = Y - Shapes->param.gen.p[1];
  816 
  817          /*  Rotate point to outer ellipse's orientation  */
  818          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  819          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  820 
  821          x /= Shapes->param.gen.p[4];
  822          y /= Shapes->param.gen.p[5];
  823          r = x*x + y*y;
  824          if( r>1.0 )
  825             comp_result = 0;
  826          else {
  827             /*  Repeat test for inner ellipse  */
  828             x =  xprime * Shapes->param.gen.b + yprime * Shapes->param.gen.a;
  829             y = -xprime * Shapes->param.gen.a + yprime * Shapes->param.gen.b;
  830 
  831             x /= Shapes->param.gen.p[2];
  832             y /= Shapes->param.gen.p[3];
  833             r = x*x + y*y;
  834             if( r<1.0 )
  835                comp_result = 0;
  836          }
  837          break;
  838 
  839       case line_rgn:
  840          /*  Shift origin to first point of line  */
  841          xprime = X - Shapes->param.gen.p[0];
  842          yprime = Y - Shapes->param.gen.p[1];
  843 
  844          /*  Rotate point to line's orientation  */
  845          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  846          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  847 
  848          if( (y < -0.5) || (y >= 0.5) || (x < -0.5)
  849              || (x >= Shapes->param.gen.a) )
  850             comp_result = 0;
  851          break;
  852 
  853       case point_rgn:
  854          /*  Shift origin to center of region  */
  855          x = X - Shapes->param.gen.p[0];
  856          y = Y - Shapes->param.gen.p[1];
  857 
  858          if ( (x<-0.5) || (x>=0.5) || (y<-0.5) || (y>=0.5) )
  859             comp_result = 0;
  860          break;
  861 
  862       case poly_rgn:
  863          if( X<Shapes->xmin || X>Shapes->xmax
  864              || Y<Shapes->ymin || Y>Shapes->ymax )
  865             comp_result = 0;
  866          else
  867             comp_result = Pt_in_Poly( X, Y, Shapes->param.poly.nPts,
  868                                        Shapes->param.poly.Pts );
  869          break;
  870 
  871       case panda_rgn:
  872          /*  Shift origin to center of region  */
  873          x = X - Shapes->param.gen.p[0];
  874          y = Y - Shapes->param.gen.p[1];
  875 
  876          r = x*x + y*y;
  877          if ( r < Shapes->param.gen.a || r > Shapes->param.gen.b ) {
  878        comp_result = 0;
  879      } else {
  880        if( x || y ) {
  881          th = atan2( y, x ) * RadToDeg;
  882          if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) {
  883                if( th < Shapes->param.gen.p[2] || th > Shapes->param.gen.p[3] )
  884          comp_result = 0;
  885          } else {
  886                if( th < Shapes->param.gen.p[2] && th > Shapes->param.gen.p[3] )
  887          comp_result = 0;
  888          }
  889        }
  890          }
  891          break;
  892 
  893       case epanda_rgn:
  894          /*  Shift origin to center of region  */
  895          xprime = X - Shapes->param.gen.p[0];
  896          yprime = Y - Shapes->param.gen.p[1];
  897 
  898          /*  Rotate point to region's orientation  */
  899          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  900          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  901      xprime = x;
  902      yprime = y;
  903 
  904      /* outer region test */
  905          x = xprime/Shapes->param.gen.p[7];
  906          y = yprime/Shapes->param.gen.p[8];
  907          r = x*x + y*y;
  908      if ( r>1.0 )
  909        comp_result = 0;
  910      else {
  911        /* inner region test */
  912        x = xprime/Shapes->param.gen.p[5];
  913        y = yprime/Shapes->param.gen.p[6];
  914        r = x*x + y*y;
  915        if ( r<1.0 )
  916          comp_result = 0;
  917        else {
  918          /* angle test */
  919          if( xprime || yprime ) {
  920            th = atan2( yprime, xprime ) * RadToDeg;
  921            if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) {
  922          if( th < Shapes->param.gen.p[2] || th > Shapes->param.gen.p[3] )
  923            comp_result = 0;
  924            } else {
  925          if( th < Shapes->param.gen.p[2] && th > Shapes->param.gen.p[3] )
  926            comp_result = 0;
  927            }
  928          }
  929        }
  930      }
  931          break;
  932 
  933       case bpanda_rgn:
  934          /*  Shift origin to center of region  */
  935          xprime = X - Shapes->param.gen.p[0];
  936          yprime = Y - Shapes->param.gen.p[1];
  937 
  938          /*  Rotate point to region's orientation  */
  939          x =  xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
  940          y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
  941 
  942      /* outer box test */
  943          dx = 0.5 * Shapes->param.gen.p[7];
  944          dy = 0.5 * Shapes->param.gen.p[8];
  945          if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
  946        comp_result = 0;
  947      else {
  948        /* inner box test */
  949        dx = 0.5 * Shapes->param.gen.p[5];
  950        dy = 0.5 * Shapes->param.gen.p[6];
  951        if( (x >= -dx) && (x <= dx) && (y >= -dy) && (y <= dy) )
  952          comp_result = 0;
  953        else {
  954          /* angle test */
  955          if( x || y ) {
  956            th = atan2( y, x ) * RadToDeg;
  957            if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) {
  958          if( th < Shapes->param.gen.p[2] || th > Shapes->param.gen.p[3] )
  959            comp_result = 0;
  960            } else {
  961          if( th < Shapes->param.gen.p[2] && th > Shapes->param.gen.p[3] )
  962            comp_result = 0;
  963            }
  964          }
  965        }
  966      }
  967          break;
  968       }
  969 
  970       if( !Shapes->sign ) comp_result = !comp_result;
  971 
  972      } 
  973 
  974    }
  975 
  976    result = result || comp_result;
  977    
  978    return( result );
  979 }
  980 
  981 /*---------------------------------------------------------------------------*/
  982 void fits_free_region( SAORegion *Rgn )
  983 /*   Free up memory allocated to hold the region data.                       
  984    This is more complicated for the case of polygons, which may be sharing
  985    points arrays due to shallow copying (in fits_set_region_components) of
  986    'exluded' regions.  We must ensure that these arrays are only freed once.       
  987 
  988 /*---------------------------------------------------------------------------*/
  989 {
  990    int i,j;
  991    
  992    int nFreedPoly=0;
  993    int nPolyArraySize=10;
  994    double **freedPolyPtrs=0;
  995    double *ptsToFree=0;
  996    int isAlreadyFreed=0;
  997    
  998    freedPolyPtrs = (double**)malloc(nPolyArraySize*sizeof(double*));
  999 
 1000    for( i=0; i<Rgn->nShapes; i++ )
 1001       if( Rgn->Shapes[i].shape == poly_rgn )
 1002       {
 1003          /* No shared arrays for 'include' polygons */
 1004          if (Rgn->Shapes[i].sign)
 1005             free(Rgn->Shapes[i].param.poly.Pts);
 1006          else
 1007          {
 1008             ptsToFree = Rgn->Shapes[i].param.poly.Pts;
 1009             isAlreadyFreed = 0;
 1010             for (j=0; j<nFreedPoly && !isAlreadyFreed; j++)
 1011             {
 1012                if (freedPolyPtrs[j] == ptsToFree)
 1013                   isAlreadyFreed = 1;
 1014             }
 1015             if (!isAlreadyFreed)
 1016             {
 1017                free(ptsToFree);
 1018                /* Now add pointer to array of freed points */
 1019                if (nFreedPoly == nPolyArraySize)
 1020                {
 1021                   nPolyArraySize *= 2;
 1022                   freedPolyPtrs = (double **)realloc(freedPolyPtrs, 
 1023                           nPolyArraySize*sizeof(double*));
 1024                }
 1025                freedPolyPtrs[nFreedPoly] = ptsToFree;
 1026                ++nFreedPoly;
 1027             }
 1028          }
 1029       }
 1030    if( Rgn->Shapes )
 1031       free( Rgn->Shapes );
 1032    free( Rgn );
 1033    
 1034    free(freedPolyPtrs);
 1035 }
 1036 
 1037 /*---------------------------------------------------------------------------*/
 1038 static int Pt_in_Poly( double x,
 1039                        double y,
 1040                        int nPts,
 1041                        double *Pts )
 1042 /*  Internal routine for testing whether the coordinate x,y is within the    */
 1043 /*  polygon region traced out by the array Pts.                              */
 1044 /*---------------------------------------------------------------------------*/
 1045 {
 1046    int i, j, flag=0;
 1047    double prevX, prevY;
 1048    double nextX, nextY;
 1049    double dx, dy, Dy;
 1050 
 1051    nextX = Pts[nPts-2];
 1052    nextY = Pts[nPts-1];
 1053 
 1054    for( i=0; i<nPts; i+=2 ) {
 1055       prevX = nextX;
 1056       prevY = nextY;
 1057 
 1058       nextX = Pts[i];
 1059       nextY = Pts[i+1];
 1060 
 1061       if( (y>prevY && y>=nextY) || (y<prevY && y<=nextY)
 1062           || (x>prevX && x>=nextX) )
 1063          continue;
 1064       
 1065       /* Check to see if x,y lies right on the segment */
 1066 
 1067       if( x>=prevX || x>nextX ) {
 1068          dy = y - prevY;
 1069          Dy = nextY - prevY;
 1070 
 1071          if( fabs(Dy)<1e-10 ) {
 1072             if( fabs(dy)<1e-10 )
 1073                return( 1 );
 1074             else
 1075                continue;
 1076          }
 1077 
 1078          dx = prevX + ( (nextX-prevX)/(Dy) ) * dy - x;
 1079          if( dx < -1e-10 )
 1080             continue;
 1081          if( dx <  1e-10 )
 1082             return( 1 );
 1083       }
 1084 
 1085       /* There is an intersection! Make sure it isn't a V point.  */
 1086 
 1087       if( y != prevY ) {
 1088          flag = 1 - flag;
 1089       } else {
 1090          j = i+1;  /* Point to Y component */
 1091          do {
 1092             if( j>1 )
 1093                j -= 2;
 1094             else
 1095                j = nPts-1;
 1096          } while( y == Pts[j] );
 1097 
 1098          if( (nextY-y)*(y-Pts[j]) > 0 )
 1099             flag = 1-flag;
 1100       }
 1101 
 1102    }
 1103    return( flag );
 1104 }
 1105 /*---------------------------------------------------------------------------*/
 1106 void fits_set_region_components ( SAORegion *aRgn )
 1107 {
 1108 /* 
 1109    Internal routine to turn a collection of regions read from an ascii file into
 1110    the more complex structure that is allowed by the FITS REGION extension with
 1111    multiple components. Regions are anded within components and ored between them
 1112    ie for a pixel to be selected it must be selected by at least one component
 1113    and to be selected by a component it must be selected by all that component's
 1114    shapes.
 1115 
 1116    The algorithm is to replicate every exclude region after every include
 1117    region before it in the list. eg reg1, reg2, -reg3, reg4, -reg5 becomes
 1118    (reg1, -reg3, -reg5), (reg2, -reg5, -reg3), (reg4, -reg5) where the
 1119    parentheses designate components.
 1120 */
 1121 
 1122   int i, j, k, icomp;
 1123 
 1124 /* loop round shapes */
 1125 
 1126   i = 0;
 1127   while ( i<aRgn->nShapes ) {
 1128 
 1129     /* first do the case of an exclude region */
 1130 
 1131     if ( !aRgn->Shapes[i].sign ) {
 1132 
 1133       /* we need to run back through the list copying the current shape as
 1134      required. start by findin the first include shape before this exclude */
 1135 
 1136       j = i-1;
 1137       while ( j > 0 && !aRgn->Shapes[j].sign ) j--;
 1138 
 1139       /* then go back one more shape */
 1140 
 1141       j--;
 1142 
 1143       /* and loop back through the regions */
 1144 
 1145       while ( j >= 0 ) {
 1146 
 1147     /* if this is an include region then insert a copy of the exclude
 1148        region immediately after it */
 1149            
 1150         /* Note that this makes shallow copies of a polygon's dynamically
 1151         allocated Pts array -- the memory is shared.  This must be checked
 1152         when freeing in fits_free_region. */
 1153 
 1154     if ( aRgn->Shapes[j].sign ) {
 1155 
 1156       aRgn->Shapes = (RgnShape *) realloc (aRgn->Shapes,(1+aRgn->nShapes)*sizeof(RgnShape));
 1157       aRgn->nShapes++;
 1158       for (k=aRgn->nShapes-1; k>j+1; k--) aRgn->Shapes[k] = aRgn->Shapes[k-1];
 1159 
 1160       i++;
 1161       aRgn->Shapes[j+1] = aRgn->Shapes[i];
 1162 
 1163     }
 1164 
 1165     j--;
 1166 
 1167       }
 1168 
 1169     }
 1170 
 1171     i++;
 1172 
 1173   }
 1174 
 1175   /* now set the component numbers */
 1176 
 1177   icomp = 0;
 1178   for ( i=0; i<aRgn->nShapes; i++ ) {
 1179     if ( aRgn->Shapes[i].sign ) icomp++;
 1180     aRgn->Shapes[i].comp = icomp;
 1181 
 1182     /*
 1183     printf("i = %d, shape = %d, sign = %d, comp = %d\n", i, aRgn->Shapes[i].shape, aRgn->Shapes[i].sign, aRgn->Shapes[i].comp);
 1184     */
 1185 
 1186   }
 1187 
 1188   return;
 1189 
 1190 }
 1191 
 1192 /*---------------------------------------------------------------------------*/
 1193 void fits_setup_shape ( RgnShape *newShape)
 1194 {
 1195 /* Perform some useful calculations now to speed up filter later             */
 1196 
 1197   double X, Y, R;
 1198   double *coords;
 1199   int i;
 1200 
 1201   if ( newShape->shape == poly_rgn ) {
 1202     coords = newShape->param.poly.Pts;
 1203   } else {
 1204     coords = newShape->param.gen.p;
 1205   }
 1206 
 1207   switch( newShape->shape ) {
 1208   case circle_rgn:
 1209     newShape->param.gen.a = coords[2] * coords[2];
 1210     break;
 1211   case annulus_rgn:
 1212     newShape->param.gen.a = coords[2] * coords[2];
 1213     newShape->param.gen.b = coords[3] * coords[3];
 1214     break;
 1215   case sector_rgn:
 1216     while( coords[2]> 180.0 ) coords[2] -= 360.0;
 1217     while( coords[2]<=-180.0 ) coords[2] += 360.0;
 1218     while( coords[3]> 180.0 ) coords[3] -= 360.0;
 1219     while( coords[3]<=-180.0 ) coords[3] += 360.0;
 1220     break;
 1221   case ellipse_rgn:
 1222     newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
 1223     newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
 1224     break;
 1225   case elliptannulus_rgn:
 1226     newShape->param.gen.a    = sin( myPI * (coords[6] / 180.0) );
 1227     newShape->param.gen.b    = cos( myPI * (coords[6] / 180.0) );
 1228     newShape->param.gen.sinT = sin( myPI * (coords[7] / 180.0) );
 1229     newShape->param.gen.cosT = cos( myPI * (coords[7] / 180.0) );
 1230     break;
 1231   case box_rgn:
 1232     newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
 1233     newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
 1234     break;
 1235   case boxannulus_rgn:
 1236     newShape->param.gen.a    = sin( myPI * (coords[6] / 180.0) );
 1237     newShape->param.gen.b    = cos( myPI * (coords[6] / 180.0) );
 1238     newShape->param.gen.sinT = sin( myPI * (coords[7] / 180.0) );
 1239     newShape->param.gen.cosT = cos( myPI * (coords[7] / 180.0) );
 1240     break;
 1241   case rectangle_rgn:
 1242     newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
 1243     newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
 1244     X = 0.5 * ( coords[2]-coords[0] );
 1245     Y = 0.5 * ( coords[3]-coords[1] );
 1246     newShape->param.gen.a = fabs( X * newShape->param.gen.cosT
 1247                   + Y * newShape->param.gen.sinT );
 1248     newShape->param.gen.b = fabs( Y * newShape->param.gen.cosT
 1249                   - X * newShape->param.gen.sinT );
 1250     newShape->param.gen.p[5] = 0.5 * ( coords[2]+coords[0] );
 1251     newShape->param.gen.p[6] = 0.5 * ( coords[3]+coords[1] );
 1252     break;
 1253   case diamond_rgn:
 1254     newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
 1255     newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
 1256     break;
 1257   case line_rgn:
 1258     X = coords[2] - coords[0];
 1259     Y = coords[3] - coords[1];
 1260     R = sqrt( X*X + Y*Y );
 1261     newShape->param.gen.sinT = ( R ? Y/R : 0.0 );
 1262     newShape->param.gen.cosT = ( R ? X/R : 1.0 );
 1263     newShape->param.gen.a    = R + 0.5;
 1264     break;
 1265   case panda_rgn:
 1266     while( coords[2]> 180.0 ) coords[2] -= 360.0;
 1267     while( coords[2]<=-180.0 ) coords[2] += 360.0;
 1268     while( coords[3]> 180.0 ) coords[3] -= 360.0;
 1269     while( coords[3]<=-180.0 ) coords[3] += 360.0;
 1270     newShape->param.gen.a = newShape->param.gen.p[5]*newShape->param.gen.p[5];
 1271     newShape->param.gen.b = newShape->param.gen.p[6]*newShape->param.gen.p[6];
 1272     break;
 1273   case epanda_rgn:
 1274   case bpanda_rgn:
 1275     while( coords[2]> 180.0 ) coords[2] -= 360.0;
 1276     while( coords[2]<=-180.0 ) coords[2] += 360.0;
 1277     while( coords[3]> 180.0 ) coords[3] -= 360.0;
 1278     while( coords[3]<=-180.0 ) coords[3] += 360.0;
 1279     newShape->param.gen.sinT = sin( myPI * (coords[10] / 180.0) );
 1280     newShape->param.gen.cosT = cos( myPI * (coords[10] / 180.0) );
 1281     break;
 1282   default:
 1283     break;
 1284   }
 1285 
 1286   /*  Set the xmin, xmax, ymin, ymax elements of the RgnShape structure */
 1287 
 1288   /* For everything which has first two parameters as center position just */
 1289   /* find a circle that encompasses the region and use it to set the       */
 1290   /* bounding box                                                          */
 1291 
 1292   R = -1.0;
 1293 
 1294   switch ( newShape->shape ) {
 1295 
 1296   case circle_rgn:
 1297     R = coords[2];
 1298     break;
 1299 
 1300   case annulus_rgn:
 1301     R = coords[3];
 1302     break;
 1303 
 1304   case ellipse_rgn:
 1305     if ( coords[2] > coords[3] ) {
 1306       R = coords[2];
 1307     } else {
 1308       R = coords[3];
 1309     }
 1310     break;
 1311 
 1312   case elliptannulus_rgn:
 1313     if ( coords[4] > coords[5] ) {
 1314       R = coords[4];
 1315     } else {
 1316       R = coords[5];
 1317     }
 1318     break;
 1319 
 1320   case box_rgn:
 1321     R = sqrt(coords[2]*coords[2]+
 1322          coords[3]*coords[3])/2.0;
 1323     break;
 1324 
 1325   case boxannulus_rgn:
 1326     R = sqrt(coords[4]*coords[5]+
 1327          coords[4]*coords[5])/2.0;
 1328     break;
 1329 
 1330   case diamond_rgn:
 1331     if ( coords[2] > coords[3] ) {
 1332       R = coords[2]/2.0;
 1333     } else {
 1334       R = coords[3]/2.0;
 1335     }
 1336     break;
 1337     
 1338   case point_rgn:
 1339     R = 1.0;
 1340     break;
 1341 
 1342   case panda_rgn:
 1343     R = coords[6];
 1344     break;
 1345 
 1346   case epanda_rgn:
 1347     if ( coords[7] > coords[8] ) {
 1348       R = coords[7];
 1349     } else {
 1350       R = coords[8];
 1351     }
 1352     break;
 1353 
 1354   case bpanda_rgn:
 1355     R = sqrt(coords[7]*coords[8]+
 1356          coords[7]*coords[8])/2.0;
 1357     break;
 1358 
 1359   default:
 1360     break;
 1361   }
 1362 
 1363   if ( R > 0.0 ) {
 1364 
 1365     newShape->xmin = coords[0] - R;
 1366     newShape->xmax = coords[0] + R;
 1367     newShape->ymin = coords[1] - R;
 1368     newShape->ymax = coords[1] + R;
 1369 
 1370     return;
 1371 
 1372   }
 1373 
 1374   /* Now do the rest of the shapes that require individual methods */
 1375 
 1376   switch ( newShape->shape ) {
 1377 
 1378   case rectangle_rgn:
 1379     R = sqrt((coords[5]-coords[0])*(coords[5]-coords[0])+
 1380          (coords[6]-coords[1])*(coords[6]-coords[1]));
 1381     newShape->xmin = coords[5] - R;
 1382     newShape->xmax = coords[5] + R;
 1383     newShape->ymin = coords[6] - R;
 1384     newShape->ymax = coords[6] + R;
 1385     break;
 1386 
 1387   case poly_rgn:
 1388     newShape->xmin = coords[0];
 1389     newShape->xmax = coords[0];
 1390     newShape->ymin = coords[1];
 1391     newShape->ymax = coords[1];
 1392     for( i=2; i < newShape->param.poly.nPts; ) {
 1393       if( newShape->xmin > coords[i] ) /* Min X */
 1394     newShape->xmin = coords[i];
 1395       if( newShape->xmax < coords[i] ) /* Max X */
 1396     newShape->xmax = coords[i];
 1397       i++;
 1398       if( newShape->ymin > coords[i] ) /* Min Y */
 1399     newShape->ymin = coords[i];
 1400       if( newShape->ymax < coords[i] ) /* Max Y */
 1401     newShape->ymax = coords[i];
 1402       i++;
 1403     }
 1404     break;
 1405 
 1406   case line_rgn:
 1407     if ( coords[0] > coords[2] ) {
 1408       newShape->xmin = coords[2];
 1409       newShape->xmax = coords[0];
 1410     } else {
 1411       newShape->xmin = coords[0];
 1412       newShape->xmax = coords[2];
 1413     }
 1414     if ( coords[1] > coords[3] ) {
 1415       newShape->ymin = coords[3];
 1416       newShape->ymax = coords[1];
 1417     } else {
 1418       newShape->ymin = coords[1];
 1419       newShape->ymax = coords[3];
 1420     }
 1421 
 1422     break;
 1423 
 1424     /* sector doesn't have min and max so indicate by setting max < min */
 1425 
 1426   case sector_rgn:
 1427     newShape->xmin = 1.0;
 1428     newShape->xmax = -1.0;
 1429     newShape->ymin = 1.0;
 1430     newShape->ymax = -1.0;
 1431     break;
 1432 
 1433   default:
 1434     break;
 1435   }
 1436 
 1437   return;
 1438 
 1439 }
 1440 
 1441 /*---------------------------------------------------------------------------*/
 1442 int fits_read_fits_region ( fitsfile *fptr, 
 1443                 WCSdata *wcs, 
 1444                 SAORegion **Rgn, 
 1445                 int *status)
 1446 /*  Read regions from a FITS region extension and return the information     */
 1447 /*  in the "SAORegion" structure.  If it is nonNULL, use wcs to convert the  */
 1448 /*  region coordinates to pixels.  Return an error if region is in degrees   */
 1449 /*  but no WCS data is provided.                                             */
 1450 /*---------------------------------------------------------------------------*/
 1451 {
 1452 
 1453   int i, j, icol[6], idum, anynul, npos;
 1454   int dotransform, got_component = 1, tstatus;
 1455   long icsize[6];
 1456   double X, Y, Theta, Xsave = 0, Ysave = 0, Xpos, Ypos;
 1457   double *coords;
 1458   char *cvalue, *cvalue2;
 1459   char comment[FLEN_COMMENT];
 1460   char colname[6][FLEN_VALUE] = {"X", "Y", "SHAPE", "R", "ROTANG", "COMPONENT"};
 1461   char shapename[17][FLEN_VALUE] = {"POINT","CIRCLE","ELLIPSE","ANNULUS",
 1462                     "ELLIPTANNULUS","BOX","ROTBOX","BOXANNULUS",
 1463                     "RECTANGLE","ROTRECTANGLE","POLYGON","PIE",
 1464                     "SECTOR","DIAMOND","RHOMBUS","ROTDIAMOND",
 1465                     "ROTRHOMBUS"};
 1466   int shapetype[17] = {point_rgn, circle_rgn, ellipse_rgn, annulus_rgn, 
 1467                elliptannulus_rgn, box_rgn, box_rgn, boxannulus_rgn, 
 1468                rectangle_rgn, rectangle_rgn, poly_rgn, sector_rgn, 
 1469                sector_rgn, diamond_rgn, diamond_rgn, diamond_rgn, 
 1470                diamond_rgn};
 1471   SAORegion *aRgn;
 1472   RgnShape *newShape;
 1473   WCSdata *regwcs = 0;
 1474 
 1475   if ( *status ) return( *status );
 1476 
 1477   aRgn = (SAORegion *)malloc( sizeof(SAORegion) );
 1478   if( ! aRgn ) {
 1479     ffpmsg("Couldn't allocate memory to hold Region file contents.");
 1480     return(*status = MEMORY_ALLOCATION );
 1481   }
 1482   aRgn->nShapes    =    0;
 1483   aRgn->Shapes     = NULL;
 1484   if( wcs && wcs->exists )
 1485     aRgn->wcs = *wcs;
 1486   else
 1487     aRgn->wcs.exists = 0;
 1488 
 1489   /* See if we are already positioned to a region extension, else */
 1490   /* move to the REGION extension (file is already open). */
 1491 
 1492   tstatus = 0;
 1493   for (i=0; i<5; i++) {
 1494     ffgcno(fptr, CASEINSEN, colname[i], &icol[i], &tstatus);
 1495   }
 1496 
 1497   if (tstatus) {
 1498     /* couldn't find the required columns, so search for "REGION" extension */
 1499     if ( ffmnhd(fptr, BINARY_TBL, "REGION", 1, status) ) {
 1500       ffpmsg("Could not move to REGION extension.");
 1501       goto error;
 1502     }
 1503   }
 1504 
 1505   /* get the number of shapes and allocate memory */
 1506 
 1507   if ( ffgky(fptr, TINT, "NAXIS2", &aRgn->nShapes, comment, status) ) {
 1508     ffpmsg("Could not read NAXIS2 keyword.");
 1509     goto error;
 1510   }
 1511 
 1512   aRgn->Shapes = (RgnShape *) malloc(aRgn->nShapes * sizeof(RgnShape));
 1513   if ( !aRgn->Shapes ) {
 1514     ffpmsg( "Failed to allocate memory for Region data");
 1515     *status = MEMORY_ALLOCATION;
 1516     goto error;
 1517   }
 1518 
 1519   /* get the required column numbers */
 1520 
 1521   for (i=0; i<5; i++) {
 1522     if ( ffgcno(fptr, CASEINSEN, colname[i], &icol[i], status) ) {
 1523       ffpmsg("Could not find column.");
 1524       goto error;
 1525     }
 1526   }
 1527 
 1528   /* try to get the optional column numbers */
 1529 
 1530   if ( ffgcno(fptr, CASEINSEN, colname[5], &icol[5], status) ) {
 1531        got_component = 0;
 1532   }
 1533 
 1534   /* if there was input WCS then read the WCS info for the region in case they */
 1535   /* are different and we have to transform */
 1536 
 1537   dotransform = 0;
 1538   if ( aRgn->wcs.exists ) {
 1539     regwcs = (WCSdata *) malloc ( sizeof(WCSdata) );
 1540     if ( !regwcs ) {
 1541       ffpmsg( "Failed to allocate memory for Region WCS data");
 1542       *status = MEMORY_ALLOCATION;
 1543       goto error;
 1544     }
 1545 
 1546     regwcs->exists = 1;
 1547     if ( ffgtcs(fptr, icol[0], icol[1], &regwcs->xrefval,  &regwcs->yrefval,
 1548         &regwcs->xrefpix, &regwcs->yrefpix, &regwcs->xinc, &regwcs->yinc,
 1549         &regwcs->rot, regwcs->type, status) ) {
 1550       regwcs->exists = 0;
 1551       *status = 0;
 1552     }
 1553 
 1554     if ( regwcs->exists && wcs->exists ) {
 1555       if ( fabs(regwcs->xrefval-wcs->xrefval) > 1.0e-6 ||
 1556        fabs(regwcs->yrefval-wcs->yrefval) > 1.0e-6 ||
 1557        fabs(regwcs->xrefpix-wcs->xrefpix) > 1.0e-6 ||
 1558        fabs(regwcs->yrefpix-wcs->yrefpix) > 1.0e-6 ||
 1559        fabs(regwcs->xinc-wcs->xinc) > 1.0e-6 ||
 1560        fabs(regwcs->yinc-wcs->yinc) > 1.0e-6 ||
 1561        fabs(regwcs->rot-wcs->rot) > 1.0e-6 ||
 1562        !strcmp(regwcs->type,wcs->type) ) dotransform = 1;
 1563     }
 1564   }
 1565 
 1566   /* get the sizes of the X, Y, R, and ROTANG vectors */
 1567 
 1568   for (i=0; i<6; i++) {
 1569     if ( ffgtdm(fptr, icol[i], 1, &idum, &icsize[i], status) ) {
 1570       ffpmsg("Could not find vector size of column.");
 1571       goto error;
 1572     }
 1573   }
 1574 
 1575   cvalue = (char *) malloc ((FLEN_VALUE+1)*sizeof(char));
 1576 
 1577   /* loop over the shapes - note 1-based counting for rows in FITS files */
 1578 
 1579   for (i=1; i<=aRgn->nShapes; i++) {
 1580 
 1581     newShape = &aRgn->Shapes[i-1];
 1582     for (j=0; j<8; j++) newShape->param.gen.p[j] = 0.0;
 1583     newShape->param.gen.a = 0.0;
 1584     newShape->param.gen.b = 0.0;
 1585     newShape->param.gen.sinT = 0.0;
 1586     newShape->param.gen.cosT = 0.0;
 1587 
 1588     /* get the shape */
 1589 
 1590     if ( ffgcvs(fptr, icol[2], i, 1, 1, " ", &cvalue, &anynul, status) ) {
 1591       ffpmsg("Could not read shape.");
 1592       goto error;
 1593     }
 1594 
 1595     /* set include or exclude */
 1596 
 1597     newShape->sign = 1;
 1598     cvalue2 = cvalue;
 1599     if ( !strncmp(cvalue,"!",1) ) {
 1600       newShape->sign = 0;
 1601       cvalue2++;
 1602     }
 1603 
 1604     /* set the shape type */
 1605 
 1606     for (j=0; j<17; j++) {
 1607       if ( !strcmp(cvalue2, shapename[j]) ) newShape->shape = shapetype[j];
 1608     }
 1609 
 1610     /* allocate memory for polygon case and set coords pointer */
 1611 
 1612     if ( newShape->shape == poly_rgn ) {
 1613       newShape->param.poly.Pts = (double *) calloc (2*icsize[0], sizeof(double));
 1614       if ( !newShape->param.poly.Pts ) {
 1615     ffpmsg("Could not allocate memory to hold polygon parameters" );
 1616     *status = MEMORY_ALLOCATION;
 1617     goto error;
 1618       }
 1619       newShape->param.poly.nPts = 2*icsize[0];
 1620       coords = newShape->param.poly.Pts;
 1621     } else {
 1622       coords = newShape->param.gen.p;
 1623     }
 1624 
 1625 
 1626   /* read X and Y. Polygon and Rectangle require special cases */
 1627 
 1628     npos = 1;
 1629     if ( newShape->shape == poly_rgn ) npos = newShape->param.poly.nPts/2;
 1630     if ( newShape->shape == rectangle_rgn ) npos = 2;
 1631 
 1632     for (j=0; j<npos; j++) {
 1633       if ( ffgcvd(fptr, icol[0], i, j+1, 1, DOUBLENULLVALUE, coords, &anynul, status) ) {
 1634     ffpmsg("Failed to read X column for polygon region");
 1635     goto error;
 1636       }
 1637       if (*coords == DOUBLENULLVALUE) {  /* check for null value end of array marker */
 1638         npos = j;
 1639     newShape->param.poly.nPts = npos * 2;
 1640     break;
 1641       }
 1642       coords++;
 1643       
 1644       if ( ffgcvd(fptr, icol[1], i, j+1, 1, DOUBLENULLVALUE, coords, &anynul, status) ) {
 1645     ffpmsg("Failed to read Y column for polygon region");
 1646     goto error;
 1647       }
 1648       if (*coords == DOUBLENULLVALUE) { /* check for null value end of array marker */
 1649         npos = j;
 1650     newShape->param.poly.nPts = npos * 2;
 1651         coords--;
 1652     break;
 1653       }
 1654       coords++;
 1655  
 1656       if (j == 0) {  /* save the first X and Y coordinate */
 1657         Xsave = *(coords - 2);
 1658     Ysave = *(coords - 1);
 1659       } else if ((Xsave == *(coords - 2)) && (Ysave == *(coords - 1)) ) {
 1660         /* if point has same coordinate as first point, this marks the end of the array */
 1661         npos = j + 1;
 1662     newShape->param.poly.nPts = npos * 2;
 1663     break;
 1664       }
 1665     }
 1666 
 1667     /* transform positions if the region and input wcs differ */
 1668 
 1669     if ( dotransform ) {
 1670 
 1671       coords -= npos*2;
 1672       Xsave = coords[0];
 1673       Ysave = coords[1];
 1674       for (j=0; j<npos; j++) {
 1675     ffwldp(coords[2*j], coords[2*j+1], regwcs->xrefval, regwcs->yrefval, regwcs->xrefpix,
 1676            regwcs->yrefpix, regwcs->xinc, regwcs->yinc, regwcs->rot,
 1677            regwcs->type, &Xpos, &Ypos, status);
 1678     ffxypx(Xpos, Ypos, wcs->xrefval, wcs->yrefval, wcs->xrefpix,
 1679            wcs->yrefpix, wcs->xinc, wcs->yinc, wcs->rot,
 1680            wcs->type, &coords[2*j], &coords[2*j+1], status);
 1681     if ( *status ) {
 1682       ffpmsg("Failed to transform coordinates");
 1683       goto error;
 1684     }
 1685       }
 1686       coords += npos*2;
 1687     }
 1688 
 1689   /* read R. Circle requires one number; Box, Diamond, Ellipse, Annulus, Sector 
 1690      and Panda two; Boxannulus and Elliptannulus four; Point, Rectangle and 
 1691      Polygon none. */
 1692 
 1693     npos = 0;
 1694     switch ( newShape->shape ) {
 1695     case circle_rgn: 
 1696       npos = 1;
 1697       break;
 1698     case box_rgn:
 1699     case diamond_rgn:
 1700     case ellipse_rgn:
 1701     case annulus_rgn:
 1702     case sector_rgn:
 1703       npos = 2;
 1704       break;
 1705     case boxannulus_rgn:
 1706     case elliptannulus_rgn:
 1707       npos = 4;
 1708       break;
 1709     default:
 1710       break;
 1711     }
 1712 
 1713     if ( npos > 0 ) {
 1714       if ( ffgcvd(fptr, icol[3], i, 1, npos, 0.0, coords, &anynul, status) ) {
 1715     ffpmsg("Failed to read R column for region");
 1716     goto error;
 1717       }
 1718 
 1719     /* transform lengths if the region and input wcs differ */
 1720 
 1721       if ( dotransform ) {
 1722     for (j=0; j<npos; j++) {
 1723       Y = Ysave + (*coords);
 1724       X = Xsave;
 1725       ffwldp(X, Y, regwcs->xrefval, regwcs->yrefval, regwcs->xrefpix,
 1726          regwcs->yrefpix, regwcs->xinc, regwcs->yinc, regwcs->rot,
 1727          regwcs->type, &Xpos, &Ypos, status);
 1728       ffxypx(Xpos, Ypos, wcs->xrefval, wcs->yrefval, wcs->xrefpix,
 1729          wcs->yrefpix, wcs->xinc, wcs->yinc, wcs->rot,
 1730          wcs->type, &X, &Y, status);
 1731       if ( *status ) {
 1732         ffpmsg("Failed to transform coordinates");
 1733         goto error;
 1734       }
 1735       *(coords++) = sqrt(pow(X-newShape->param.gen.p[0],2)+pow(Y-newShape->param.gen.p[1],2));
 1736     }
 1737       } else {
 1738     coords += npos;
 1739       }
 1740     }
 1741 
 1742   /* read ROTANG. Requires two values for Boxannulus, Elliptannulus, Sector, 
 1743      Panda; one for Box, Diamond, Ellipse; and none for Circle, Point, Annulus, 
 1744      Rectangle, Polygon */
 1745 
 1746     npos = 0;
 1747     switch ( newShape->shape ) {
 1748     case box_rgn:
 1749     case diamond_rgn:
 1750     case ellipse_rgn:
 1751       npos = 1;
 1752       break;
 1753     case boxannulus_rgn:
 1754     case elliptannulus_rgn:
 1755     case sector_rgn:
 1756       npos = 2;
 1757       break;
 1758     default:
 1759      break;
 1760     }
 1761 
 1762     if ( npos > 0 ) {
 1763       if ( ffgcvd(fptr, icol[4], i, 1, npos, 0.0, coords, &anynul, status) ) {
 1764     ffpmsg("Failed to read ROTANG column for region");
 1765     goto error;
 1766       }
 1767 
 1768     /* transform angles if the region and input wcs differ */
 1769 
 1770       if ( dotransform ) {
 1771     Theta = (wcs->rot) - (regwcs->rot);
 1772     for (j=0; j<npos; j++) *(coords++) += Theta;
 1773       } else {
 1774     coords += npos;
 1775       }
 1776     }
 1777 
 1778   /* read the component number */
 1779 
 1780     if (got_component) {
 1781       if ( ffgcv(fptr, TINT, icol[5], i, 1, 1, 0, &newShape->comp, &anynul, status) ) {
 1782         ffpmsg("Failed to read COMPONENT column for region");
 1783         goto error;
 1784       }
 1785     } else {
 1786       newShape->comp = 1;
 1787     }
 1788 
 1789 
 1790     /* do some precalculations to speed up tests */
 1791 
 1792     fits_setup_shape(newShape);
 1793 
 1794     /* end loop over shapes */
 1795 
 1796   }
 1797 
 1798 error:
 1799 
 1800    if( *status )
 1801       fits_free_region( aRgn );
 1802    else
 1803       *Rgn = aRgn;
 1804 
 1805    ffclos(fptr, status);
 1806 
 1807    return( *status );
 1808 }
 1809