From fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 Mon Sep 17 00:00:00 2001 From: Joseph Hunkeler Date: Wed, 8 Jul 2015 20:46:52 -0400 Subject: Initial commit --- vendor/cfitsio/region.c | 1747 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 1747 insertions(+) create mode 100644 vendor/cfitsio/region.c (limited to 'vendor/cfitsio/region.c') diff --git a/vendor/cfitsio/region.c b/vendor/cfitsio/region.c new file mode 100644 index 00000000..3ec5bc28 --- /dev/null +++ b/vendor/cfitsio/region.c @@ -0,0 +1,1747 @@ +#include +#include +#include +#include +#include +#include "fitsio2.h" +#include "region.h" +static int Pt_in_Poly( double x, double y, int nPts, double *Pts ); + +/*---------------------------------------------------------------------------*/ +int fits_read_rgnfile( const char *filename, + WCSdata *wcs, + SAORegion **Rgn, + int *status ) +/* Read regions from either a FITS or ASCII 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. */ +/*---------------------------------------------------------------------------*/ +{ + fitsfile *fptr; + int tstatus = 0; + + if( *status ) return( *status ); + + /* try to open as a FITS file - if that doesn't work treat as an ASCII file */ + + fits_write_errmark(); + if ( ffopen(&fptr, filename, READONLY, &tstatus) ) { + fits_clear_errmark(); + fits_read_ascii_region(filename, wcs, Rgn, status); + } else { + fits_read_fits_region(fptr, wcs, Rgn, status); + } + + return(*status); + +} +/*---------------------------------------------------------------------------*/ +int fits_read_ascii_region( 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, X, Y, x, y, ss, div, 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( isspace(*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( isspace(*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++; + if ( paramPtr ) cFmt = hhmmss_fmt; /* set format if parameter has : */ + break; + case 'd': + currLoc++; + if ( paramPtr ) cFmt = degree_fmt; /* set format if parameter has d */ + 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( isspace(*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, "linear;", 7 ) ) { + namePtr += 7; + 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 5 cases support region files created by POW + (or ds9 Version 4.x) 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, "physical", 8 ) ) { + 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; + for (i=0; i<8; i++) newShape->param.gen.p[i] = 0.0; + newShape->param.gen.a = 0.0; + newShape->param.gen.b = 0.0; + newShape->param.gen.sinT = 0.0; + newShape->param.gen.cosT = 0.0; + + while( isspace(*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( isspace(*namePtr) ) namePtr++; + if( *namePtr=='\0' ) { + ffpmsg( "Error reading Region file" ); + *status = PARSE_SYNTAX_ERR; + goto error; + } + lineLen = strlen( namePtr ) - 1; + while( isspace(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" ) ) { + if( nParams < 4 || nParams > 8 ) { + *status = PARSE_SYNTAX_ERR; + } else if ( nParams < 6 ) { + newShape->shape = ellipse_rgn; + newShape->param.gen.p[4] = 0.0; + } else { + newShape->shape = elliptannulus_rgn; + newShape->param.gen.p[6] = 0.0; + newShape->param.gen.p[7] = 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" ) ) { + if( nParams < 4 || nParams > 8 ) { + *status = PARSE_SYNTAX_ERR; + } else if ( nParams < 6 ) { + newShape->shape = box_rgn; + newShape->param.gen.p[4] = 0.0; + } else { + newShape->shape = boxannulus_rgn; + newShape->param.gen.p[6] = 0.0; + newShape->param.gen.p[7] = 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 if( !strcasecmp( namePtr, "panda" ) ) { + newShape->shape = panda_rgn; + if( nParams != 8 ) + *status = PARSE_SYNTAX_ERR; + nCoords = 2; + } else if( !strcasecmp( namePtr, "epanda" ) ) { + newShape->shape = epanda_rgn; + if( nParams < 10 || nParams > 11 ) + *status = PARSE_SYNTAX_ERR; + newShape->param.gen.p[10] = 0.0; + nCoords = 2; + } else if( !strcasecmp( namePtr, "bpanda" ) ) { + newShape->shape = bpanda_rgn; + if( nParams < 10 || nParams > 11 ) + *status = PARSE_SYNTAX_ERR; + newShape->param.gen.p[10] = 0.0; + nCoords = 2; + } 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; iexists ) { + 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( ; ixrefval, 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) ); + + } + } + + /* special case for elliptannulus and boxannulus if only one angle + was given */ + + if ( (newShape->shape == elliptannulus_rgn || + newShape->shape == boxannulus_rgn ) && nParams == 7 ) { + coords[7] = coords[6]; + } + + /* Also, correct the position angle for any WCS rotation: */ + /* If regions are specified in WCS coordintes, then the angles */ + /* are relative to the WCS system, not the pixel X,Y system */ + + if( cFmt!=pixel_fmt ) { + switch( newShape->shape ) { + case sector_rgn: + case panda_rgn: + coords[2] += (wcs->rot); + coords[3] += (wcs->rot); + break; + case box_rgn: + case rectangle_rgn: + case diamond_rgn: + case ellipse_rgn: + coords[4] += (wcs->rot); + break; + case boxannulus_rgn: + case elliptannulus_rgn: + coords[6] += (wcs->rot); + coords[7] += (wcs->rot); + break; + case epanda_rgn: + case bpanda_rgn: + coords[2] += (wcs->rot); + coords[3] += (wcs->rot); + coords[10] += (wcs->rot); + } + } + + /* do some precalculations to speed up tests */ + + fits_setup_shape(newShape); + + } /* 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) ) */ + + /* set up component numbers */ + + fits_set_region_components( aRgn ); + +error: + + if( *status ) { + fits_free_region( aRgn ); + } else { + *Rgn = aRgn; + } + + fclose( rgnFile ); + free( currLine ); + + return( *status ); +} + +/*---------------------------------------------------------------------------*/ +int fits_in_region( 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, th; + RgnShape *Shapes; + int i, cur_comp; + int result, comp_result; + + Shapes = Rgn->Shapes; + + result = 0; + comp_result = 0; + cur_comp = Rgn->Shapes[0].comp; + + for( i=0; inShapes; i++, Shapes++ ) { + + /* if this region has a different component number to the last one */ + /* then replace the accumulated selection logical with the union of */ + /* the current logical and the total logical. Reinitialize the */ + /* temporary logical. */ + + if ( i==0 || Shapes->comp != cur_comp ) { + result = result || comp_result; + cur_comp = Shapes->comp; + /* if an excluded region is given first, then implicitly */ + /* assume a previous shape that includes the entire image. */ + comp_result = !Shapes->sign; + } + + /* 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 ( (!comp_result && Shapes->sign) || (comp_result && !Shapes->sign) ) { + + comp_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) ) + comp_result = 0; + break; + + case boxannulus_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[4]; + dy = 0.5 * Shapes->param.gen.p[5]; + if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) ) { + comp_result = 0; + } else { + /* Repeat test for inner box */ + x = xprime * Shapes->param.gen.b + yprime * Shapes->param.gen.a; + y = -xprime * Shapes->param.gen.a + yprime * Shapes->param.gen.b; + + 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) ) + comp_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) ) + comp_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 ) + comp_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 ) + comp_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 ) + comp_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 ) * RadToDeg; + if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) { + if( r < Shapes->param.gen.p[2] || r > Shapes->param.gen.p[3] ) + comp_result = 0; + } else { + if( r < Shapes->param.gen.p[2] && r > Shapes->param.gen.p[3] ) + comp_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 ) + comp_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 ) + comp_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 ) + comp_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) ) + comp_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) ) + comp_result = 0; + break; + + case poly_rgn: + if( Xxmin || X>Shapes->xmax + || Yymin || Y>Shapes->ymax ) + comp_result = 0; + else + comp_result = Pt_in_Poly( X, Y, Shapes->param.poly.nPts, + Shapes->param.poly.Pts ); + break; + + case panda_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 ) { + comp_result = 0; + } else { + if( x || y ) { + th = atan2( y, x ) * RadToDeg; + if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) { + if( th < Shapes->param.gen.p[2] || th > Shapes->param.gen.p[3] ) + comp_result = 0; + } else { + if( th < Shapes->param.gen.p[2] && th > Shapes->param.gen.p[3] ) + comp_result = 0; + } + } + } + break; + + case epanda_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; + xprime = x; + yprime = y; + + /* outer region test */ + x = xprime/Shapes->param.gen.p[7]; + y = yprime/Shapes->param.gen.p[8]; + r = x*x + y*y; + if ( r>1.0 ) + comp_result = 0; + else { + /* inner region test */ + x = xprime/Shapes->param.gen.p[5]; + y = yprime/Shapes->param.gen.p[6]; + r = x*x + y*y; + if ( r<1.0 ) + comp_result = 0; + else { + /* angle test */ + if( xprime || yprime ) { + th = atan2( yprime, xprime ) * RadToDeg; + if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) { + if( th < Shapes->param.gen.p[2] || th > Shapes->param.gen.p[3] ) + comp_result = 0; + } else { + if( th < Shapes->param.gen.p[2] && th > Shapes->param.gen.p[3] ) + comp_result = 0; + } + } + } + } + break; + + case bpanda_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; + + /* outer box test */ + dx = 0.5 * Shapes->param.gen.p[7]; + dy = 0.5 * Shapes->param.gen.p[8]; + if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) ) + comp_result = 0; + else { + /* inner box test */ + dx = 0.5 * Shapes->param.gen.p[5]; + dy = 0.5 * Shapes->param.gen.p[6]; + if( (x >= -dx) && (x <= dx) && (y >= -dy) && (y <= dy) ) + comp_result = 0; + else { + /* angle test */ + if( x || y ) { + th = atan2( y, x ) * RadToDeg; + if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) { + if( th < Shapes->param.gen.p[2] || th > Shapes->param.gen.p[3] ) + comp_result = 0; + } else { + if( th < Shapes->param.gen.p[2] && th > Shapes->param.gen.p[3] ) + comp_result = 0; + } + } + } + } + break; + } + + if( !Shapes->sign ) comp_result = !comp_result; + + } + + } + + result = result || comp_result; + + return( result ); +} + +/*---------------------------------------------------------------------------*/ +void fits_free_region( SAORegion *Rgn ) +/* Free up memory allocated to hold the region data. */ +/*---------------------------------------------------------------------------*/ +{ + int i; + + for( i=0; inShapes; 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; iprevY && y>=nextY) || (yprevX && 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 ); +} +/*---------------------------------------------------------------------------*/ +void fits_set_region_components ( SAORegion *aRgn ) +{ +/* + Internal routine to turn a collection of regions read from an ascii file into + the more complex structure that is allowed by the FITS REGION extension with + multiple components. Regions are anded within components and ored between them + ie for a pixel to be selected it must be selected by at least one component + and to be selected by a component it must be selected by all that component's + shapes. + + The algorithm is to replicate every exclude region after every include + region before it in the list. eg reg1, reg2, -reg3, reg4, -reg5 becomes + (reg1, -reg3, -reg5), (reg2, -reg5, -reg3), (reg4, -reg5) where the + parentheses designate components. +*/ + + int i, j, k, icomp; + +/* loop round shapes */ + + i = 0; + while ( inShapes ) { + + /* first do the case of an exclude region */ + + if ( !aRgn->Shapes[i].sign ) { + + /* we need to run back through the list copying the current shape as + required. start by findin the first include shape before this exclude */ + + j = i-1; + while ( j > 0 && !aRgn->Shapes[j].sign ) j--; + + /* then go back one more shape */ + + j--; + + /* and loop back through the regions */ + + while ( j >= 0 ) { + + /* if this is an include region then insert a copy of the exclude + region immediately after it */ + + if ( aRgn->Shapes[j].sign ) { + + aRgn->Shapes = (RgnShape *) realloc (aRgn->Shapes,(1+aRgn->nShapes)*sizeof(RgnShape)); + aRgn->nShapes++; + for (k=aRgn->nShapes-1; k>j+1; k--) aRgn->Shapes[k] = aRgn->Shapes[k-1]; + + i++; + aRgn->Shapes[j+1] = aRgn->Shapes[i]; + + } + + j--; + + } + + } + + i++; + + } + + /* now set the component numbers */ + + icomp = 0; + for ( i=0; inShapes; i++ ) { + if ( aRgn->Shapes[i].sign ) icomp++; + aRgn->Shapes[i].comp = icomp; + + /* + printf("i = %d, shape = %d, sign = %d, comp = %d\n", i, aRgn->Shapes[i].shape, aRgn->Shapes[i].sign, aRgn->Shapes[i].comp); + */ + + } + + return; + +} + +/*---------------------------------------------------------------------------*/ +void fits_setup_shape ( RgnShape *newShape) +{ +/* Perform some useful calculations now to speed up filter later */ + + double X, Y, R; + double *coords; + int i; + + if ( newShape->shape == poly_rgn ) { + coords = newShape->param.poly.Pts; + } else { + coords = newShape->param.gen.p; + } + + 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 boxannulus_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 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 panda_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; + newShape->param.gen.a = newShape->param.gen.p[5]*newShape->param.gen.p[5]; + newShape->param.gen.b = newShape->param.gen.p[6]*newShape->param.gen.p[6]; + break; + case epanda_rgn: + case bpanda_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; + newShape->param.gen.sinT = sin( myPI * (coords[10] / 180.0) ); + newShape->param.gen.cosT = cos( myPI * (coords[10] / 180.0) ); + break; + } + + /* Set the xmin, xmax, ymin, ymax elements of the RgnShape structure */ + + /* For everything which has first two parameters as center position just */ + /* find a circle that encompasses the region and use it to set the */ + /* bounding box */ + + R = -1.0; + + switch ( newShape->shape ) { + + case circle_rgn: + R = coords[2]; + break; + + case annulus_rgn: + R = coords[3]; + break; + + case ellipse_rgn: + if ( coords[2] > coords[3] ) { + R = coords[2]; + } else { + R = coords[3]; + } + break; + + case elliptannulus_rgn: + if ( coords[4] > coords[5] ) { + R = coords[4]; + } else { + R = coords[5]; + } + break; + + case box_rgn: + R = sqrt(coords[2]*coords[2]+ + coords[3]*coords[3])/2.0; + break; + + case boxannulus_rgn: + R = sqrt(coords[4]*coords[5]+ + coords[4]*coords[5])/2.0; + break; + + case diamond_rgn: + if ( coords[2] > coords[3] ) { + R = coords[2]/2.0; + } else { + R = coords[3]/2.0; + } + break; + + case point_rgn: + R = 1.0; + break; + + case panda_rgn: + R = coords[6]; + break; + + case epanda_rgn: + if ( coords[7] > coords[8] ) { + R = coords[7]; + } else { + R = coords[8]; + } + break; + + case bpanda_rgn: + R = sqrt(coords[7]*coords[8]+ + coords[7]*coords[8])/2.0; + break; + + } + + if ( R > 0.0 ) { + + newShape->xmin = coords[0] - R; + newShape->xmax = coords[0] + R; + newShape->ymin = coords[1] - R; + newShape->ymax = coords[1] + R; + + return; + + } + + /* Now do the rest of the shapes that require individual methods */ + + switch ( newShape->shape ) { + + case rectangle_rgn: + R = sqrt((coords[5]-coords[0])*(coords[5]-coords[0])+ + (coords[6]-coords[1])*(coords[6]-coords[1])); + newShape->xmin = coords[5] - R; + newShape->xmax = coords[5] + R; + newShape->ymin = coords[6] - R; + newShape->ymax = coords[6] + R; + break; + + case poly_rgn: + newShape->xmin = coords[0]; + newShape->xmax = coords[0]; + newShape->ymin = coords[1]; + newShape->ymax = coords[1]; + for( i=2; i < newShape->param.poly.nPts; ) { + if( newShape->xmin > coords[i] ) /* Min X */ + newShape->xmin = coords[i]; + if( newShape->xmax < coords[i] ) /* Max X */ + newShape->xmax = coords[i]; + i++; + if( newShape->ymin > coords[i] ) /* Min Y */ + newShape->ymin = coords[i]; + if( newShape->ymax < coords[i] ) /* Max Y */ + newShape->ymax = coords[i]; + i++; + } + break; + + case line_rgn: + if ( coords[0] > coords[2] ) { + newShape->xmin = coords[2]; + newShape->xmax = coords[0]; + } else { + newShape->xmin = coords[0]; + newShape->xmax = coords[2]; + } + if ( coords[1] > coords[3] ) { + newShape->ymin = coords[3]; + newShape->ymax = coords[1]; + } else { + newShape->ymin = coords[1]; + newShape->ymax = coords[3]; + } + + break; + + /* sector doesn't have min and max so indicate by setting max < min */ + + case sector_rgn: + newShape->xmin = 1.0; + newShape->xmax = -1.0; + newShape->ymin = 1.0; + newShape->ymax = -1.0; + break; + + } + + return; + +} + +/*---------------------------------------------------------------------------*/ +int fits_read_fits_region ( fitsfile *fptr, + WCSdata *wcs, + SAORegion **Rgn, + int *status) +/* Read regions from a FITS region extension 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. */ +/*---------------------------------------------------------------------------*/ +{ + + int i, j, icol[6], idum, anynul, npos; + int dotransform, got_component = 1, tstatus; + long icsize[6]; + double X, Y, Theta, Xsave, Ysave, Xpos, Ypos; + double *coords; + char *cvalue, *cvalue2; + char comment[FLEN_COMMENT]; + char colname[6][FLEN_VALUE] = {"X", "Y", "SHAPE", "R", "ROTANG", "COMPONENT"}; + char shapename[17][FLEN_VALUE] = {"POINT","CIRCLE","ELLIPSE","ANNULUS", + "ELLIPTANNULUS","BOX","ROTBOX","BOXANNULUS", + "RECTANGLE","ROTRECTANGLE","POLYGON","PIE", + "SECTOR","DIAMOND","RHOMBUS","ROTDIAMOND", + "ROTRHOMBUS"}; + int shapetype[17] = {point_rgn, circle_rgn, ellipse_rgn, annulus_rgn, + elliptannulus_rgn, box_rgn, box_rgn, boxannulus_rgn, + rectangle_rgn, rectangle_rgn, poly_rgn, sector_rgn, + sector_rgn, diamond_rgn, diamond_rgn, diamond_rgn, + diamond_rgn}; + SAORegion *aRgn; + RgnShape *newShape; + WCSdata *regwcs; + + 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; + + /* See if we are already positioned to a region extension, else */ + /* move to the REGION extension (file is already open). */ + + tstatus = 0; + for (i=0; i<5; i++) { + ffgcno(fptr, CASEINSEN, colname[i], &icol[i], &tstatus); + } + + if (tstatus) { + /* couldn't find the required columns, so search for "REGION" extension */ + if ( ffmnhd(fptr, BINARY_TBL, "REGION", 1, status) ) { + ffpmsg("Could not move to REGION extension."); + goto error; + } + } + + /* get the number of shapes and allocate memory */ + + if ( ffgky(fptr, TINT, "NAXIS2", &aRgn->nShapes, comment, status) ) { + ffpmsg("Could not read NAXIS2 keyword."); + goto error; + } + + aRgn->Shapes = (RgnShape *) malloc(aRgn->nShapes * sizeof(RgnShape)); + if ( !aRgn->Shapes ) { + ffpmsg( "Failed to allocate memory for Region data"); + *status = MEMORY_ALLOCATION; + goto error; + } + + /* get the required column numbers */ + + for (i=0; i<5; i++) { + if ( ffgcno(fptr, CASEINSEN, colname[i], &icol[i], status) ) { + ffpmsg("Could not find column."); + goto error; + } + } + + /* try to get the optional column numbers */ + + if ( ffgcno(fptr, CASEINSEN, colname[5], &icol[5], status) ) { + got_component = 0; + } + + /* if there was input WCS then read the WCS info for the region in case they */ + /* are different and we have to transform */ + + dotransform = 0; + if ( aRgn->wcs.exists ) { + regwcs = (WCSdata *) malloc ( sizeof(WCSdata) ); + if ( !regwcs ) { + ffpmsg( "Failed to allocate memory for Region WCS data"); + *status = MEMORY_ALLOCATION; + goto error; + } + + regwcs->exists = 1; + if ( ffgtcs(fptr, icol[0], icol[1], ®wcs->xrefval, ®wcs->yrefval, + ®wcs->xrefpix, ®wcs->yrefpix, ®wcs->xinc, ®wcs->yinc, + ®wcs->rot, regwcs->type, status) ) { + regwcs->exists = 0; + *status = 0; + } + + if ( regwcs->exists && wcs->exists ) { + if ( fabs(regwcs->xrefval-wcs->xrefval) > 1.0e-6 || + fabs(regwcs->yrefval-wcs->yrefval) > 1.0e-6 || + fabs(regwcs->xrefpix-wcs->xrefpix) > 1.0e-6 || + fabs(regwcs->yrefpix-wcs->yrefpix) > 1.0e-6 || + fabs(regwcs->xinc-wcs->xinc) > 1.0e-6 || + fabs(regwcs->yinc-wcs->yinc) > 1.0e-6 || + fabs(regwcs->rot-wcs->rot) > 1.0e-6 || + !strcmp(regwcs->type,wcs->type) ) dotransform = 1; + } + } + + /* get the sizes of the X, Y, R, and ROTANG vectors */ + + for (i=0; i<6; i++) { + if ( ffgtdm(fptr, icol[i], 1, &idum, &icsize[i], status) ) { + ffpmsg("Could not find vector size of column."); + goto error; + } + } + + cvalue = (char *) malloc ((FLEN_VALUE+1)*sizeof(char)); + + /* loop over the shapes - note 1-based counting for rows in FITS files */ + + for (i=1; i<=aRgn->nShapes; i++) { + + newShape = &aRgn->Shapes[i-1]; + for (j=0; j<8; j++) newShape->param.gen.p[j] = 0.0; + newShape->param.gen.a = 0.0; + newShape->param.gen.b = 0.0; + newShape->param.gen.sinT = 0.0; + newShape->param.gen.cosT = 0.0; + + /* get the shape */ + + if ( ffgcvs(fptr, icol[2], i, 1, 1, " ", &cvalue, &anynul, status) ) { + ffpmsg("Could not read shape."); + goto error; + } + + /* set include or exclude */ + + newShape->sign = 1; + cvalue2 = cvalue; + if ( !strncmp(cvalue,"!",1) ) { + newShape->sign = 0; + cvalue2++; + } + + /* set the shape type */ + + for (j=0; j<9; j++) { + if ( !strcmp(cvalue2, shapename[j]) ) newShape->shape = shapetype[j]; + } + + /* allocate memory for polygon case and set coords pointer */ + + if ( newShape->shape == poly_rgn ) { + newShape->param.poly.Pts = (double *) calloc (2*icsize[0], 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 = 2*icsize[0]; + coords = newShape->param.poly.Pts; + } else { + coords = newShape->param.gen.p; + } + + + /* read X and Y. Polygon and Rectangle require special cases */ + + npos = 1; + if ( newShape->shape == poly_rgn ) npos = newShape->param.poly.nPts/2; + if ( newShape->shape == rectangle_rgn ) npos = 2; + + for (j=0; jparam.poly.nPts = npos * 2; + break; + } + coords++; + + if ( ffgcvd(fptr, icol[1], i, j+1, 1, DOUBLENULLVALUE, coords, &anynul, status) ) { + ffpmsg("Failed to read Y column for polygon region"); + goto error; + } + if (*coords == DOUBLENULLVALUE) { /* check for null value end of array marker */ + npos = j; + newShape->param.poly.nPts = npos * 2; + coords--; + break; + } + coords++; + + if (j == 0) { /* save the first X and Y coordinate */ + Xsave = *(coords - 2); + Ysave = *(coords - 1); + } else if ((Xsave == *(coords - 2)) && (Ysave == *(coords - 1)) ) { + /* if point has same coordinate as first point, this marks the end of the array */ + npos = j + 1; + newShape->param.poly.nPts = npos * 2; + break; + } + } + + /* transform positions if the region and input wcs differ */ + + if ( dotransform ) { + + coords -= npos*2; + Xsave = coords[0]; + Ysave = coords[1]; + for (j=0; jxrefval, regwcs->yrefval, regwcs->xrefpix, + regwcs->yrefpix, regwcs->xinc, regwcs->yinc, regwcs->rot, + regwcs->type, &Xpos, &Ypos, status); + ffxypx(Xpos, Ypos, wcs->xrefval, wcs->yrefval, wcs->xrefpix, + wcs->yrefpix, wcs->xinc, wcs->yinc, wcs->rot, + wcs->type, &coords[2*j], &coords[2*j+1], status); + if ( *status ) { + ffpmsg("Failed to transform coordinates"); + goto error; + } + } + coords += npos*2; + } + + /* read R. Circle requires one number; Box, Diamond, Ellipse, Annulus, Sector + and Panda two; Boxannulus and Elliptannulus four; Point, Rectangle and + Polygon none. */ + + npos = 0; + switch ( newShape->shape ) { + case circle_rgn: + npos = 1; + break; + case box_rgn: + case diamond_rgn: + case ellipse_rgn: + case annulus_rgn: + case sector_rgn: + npos = 2; + break; + case boxannulus_rgn: + case elliptannulus_rgn: + npos = 4; + break; + } + + if ( npos > 0 ) { + if ( ffgcvd(fptr, icol[3], i, 1, npos, 0.0, coords, &anynul, status) ) { + ffpmsg("Failed to read R column for region"); + goto error; + } + + /* transform lengths if the region and input wcs differ */ + + if ( dotransform ) { + for (j=0; jxrefval, regwcs->yrefval, regwcs->xrefpix, + regwcs->yrefpix, regwcs->xinc, regwcs->yinc, regwcs->rot, + regwcs->type, &Xpos, &Ypos, status); + ffxypx(Xpos, Ypos, wcs->xrefval, wcs->yrefval, wcs->xrefpix, + wcs->yrefpix, wcs->xinc, wcs->yinc, wcs->rot, + wcs->type, &X, &Y, status); + if ( *status ) { + ffpmsg("Failed to transform coordinates"); + goto error; + } + *(coords++) = sqrt(pow(X-newShape->param.gen.p[0],2)+pow(Y-newShape->param.gen.p[1],2)); + } + } else { + coords += npos; + } + } + + /* read ROTANG. Requires two values for Boxannulus, Elliptannulus, Sector, + Panda; one for Box, Diamond, Ellipse; and none for Circle, Point, Annulus, + Rectangle, Polygon */ + + npos = 0; + switch ( newShape->shape ) { + case box_rgn: + case diamond_rgn: + case ellipse_rgn: + npos = 1; + break; + case boxannulus_rgn: + case elliptannulus_rgn: + case sector_rgn: + npos = 2; + break; + } + + if ( npos > 0 ) { + if ( ffgcvd(fptr, icol[4], i, 1, npos, 0.0, coords, &anynul, status) ) { + ffpmsg("Failed to read ROTANG column for region"); + goto error; + } + + /* transform angles if the region and input wcs differ */ + + if ( dotransform ) { + Theta = (wcs->rot) - (regwcs->rot); + for (j=0; jcomp, &anynul, status) ) { + ffpmsg("Failed to read COMPONENT column for region"); + goto error; + } + } else { + newShape->comp = 1; + } + + + /* do some precalculations to speed up tests */ + + fits_setup_shape(newShape); + + /* end loop over shapes */ + + } + +error: + + if( *status ) + fits_free_region( aRgn ); + else + *Rgn = aRgn; + + ffclos(fptr, status); + + return( *status ); +} + -- cgit