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 /*---------------------------------------------------------------------------*/
fits_read_rgnfile(const char * filename,WCSdata * wcs,SAORegion ** Rgn,int * status)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 /*---------------------------------------------------------------------------*/
fits_read_ascii_region(const char * filename,WCSdata * wcs,SAORegion ** Rgn,int * status)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 /*---------------------------------------------------------------------------*/
fits_in_region(double X,double Y,SAORegion * Rgn)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 /*---------------------------------------------------------------------------*/
fits_free_region(SAORegion * Rgn)982 void fits_free_region( SAORegion *Rgn )
983 /* Free up memory allocated to hold the region data. */
984 /*---------------------------------------------------------------------------*/
985 {
986 int i;
987
988 for( i=0; i<Rgn->nShapes; i++ )
989 if( Rgn->Shapes[i].shape == poly_rgn )
990 free( Rgn->Shapes[i].param.poly.Pts );
991 if( Rgn->Shapes )
992 free( Rgn->Shapes );
993 free( Rgn );
994 }
995
996 /*---------------------------------------------------------------------------*/
Pt_in_Poly(double x,double y,int nPts,double * Pts)997 static int Pt_in_Poly( double x,
998 double y,
999 int nPts,
1000 double *Pts )
1001 /* Internal routine for testing whether the coordinate x,y is within the */
1002 /* polygon region traced out by the array Pts. */
1003 /*---------------------------------------------------------------------------*/
1004 {
1005 int i, j, flag=0;
1006 double prevX, prevY;
1007 double nextX, nextY;
1008 double dx, dy, Dy;
1009
1010 nextX = Pts[nPts-2];
1011 nextY = Pts[nPts-1];
1012
1013 for( i=0; i<nPts; i+=2 ) {
1014 prevX = nextX;
1015 prevY = nextY;
1016
1017 nextX = Pts[i];
1018 nextY = Pts[i+1];
1019
1020 if( (y>prevY && y>=nextY) || (y<prevY && y<=nextY)
1021 || (x>prevX && x>=nextX) )
1022 continue;
1023
1024 /* Check to see if x,y lies right on the segment */
1025
1026 if( x>=prevX || x>nextX ) {
1027 dy = y - prevY;
1028 Dy = nextY - prevY;
1029
1030 if( fabs(Dy)<1e-10 ) {
1031 if( fabs(dy)<1e-10 )
1032 return( 1 );
1033 else
1034 continue;
1035 }
1036
1037 dx = prevX + ( (nextX-prevX)/(Dy) ) * dy - x;
1038 if( dx < -1e-10 )
1039 continue;
1040 if( dx < 1e-10 )
1041 return( 1 );
1042 }
1043
1044 /* There is an intersection! Make sure it isn't a V point. */
1045
1046 if( y != prevY ) {
1047 flag = 1 - flag;
1048 } else {
1049 j = i+1; /* Point to Y component */
1050 do {
1051 if( j>1 )
1052 j -= 2;
1053 else
1054 j = nPts-1;
1055 } while( y == Pts[j] );
1056
1057 if( (nextY-y)*(y-Pts[j]) > 0 )
1058 flag = 1-flag;
1059 }
1060
1061 }
1062 return( flag );
1063 }
1064 /*---------------------------------------------------------------------------*/
fits_set_region_components(SAORegion * aRgn)1065 void fits_set_region_components ( SAORegion *aRgn )
1066 {
1067 /*
1068 Internal routine to turn a collection of regions read from an ascii file into
1069 the more complex structure that is allowed by the FITS REGION extension with
1070 multiple components. Regions are anded within components and ored between them
1071 ie for a pixel to be selected it must be selected by at least one component
1072 and to be selected by a component it must be selected by all that component's
1073 shapes.
1074
1075 The algorithm is to replicate every exclude region after every include
1076 region before it in the list. eg reg1, reg2, -reg3, reg4, -reg5 becomes
1077 (reg1, -reg3, -reg5), (reg2, -reg5, -reg3), (reg4, -reg5) where the
1078 parentheses designate components.
1079 */
1080
1081 int i, j, k, icomp;
1082
1083 /* loop round shapes */
1084
1085 i = 0;
1086 while ( i<aRgn->nShapes ) {
1087
1088 /* first do the case of an exclude region */
1089
1090 if ( !aRgn->Shapes[i].sign ) {
1091
1092 /* we need to run back through the list copying the current shape as
1093 required. start by findin the first include shape before this exclude */
1094
1095 j = i-1;
1096 while ( j > 0 && !aRgn->Shapes[j].sign ) j--;
1097
1098 /* then go back one more shape */
1099
1100 j--;
1101
1102 /* and loop back through the regions */
1103
1104 while ( j >= 0 ) {
1105
1106 /* if this is an include region then insert a copy of the exclude
1107 region immediately after it */
1108
1109 if ( aRgn->Shapes[j].sign ) {
1110
1111 aRgn->Shapes = (RgnShape *) realloc (aRgn->Shapes,(1+aRgn->nShapes)*sizeof(RgnShape));
1112 aRgn->nShapes++;
1113 for (k=aRgn->nShapes-1; k>j+1; k--) aRgn->Shapes[k] = aRgn->Shapes[k-1];
1114
1115 i++;
1116 aRgn->Shapes[j+1] = aRgn->Shapes[i];
1117
1118 }
1119
1120 j--;
1121
1122 }
1123
1124 }
1125
1126 i++;
1127
1128 }
1129
1130 /* now set the component numbers */
1131
1132 icomp = 0;
1133 for ( i=0; i<aRgn->nShapes; i++ ) {
1134 if ( aRgn->Shapes[i].sign ) icomp++;
1135 aRgn->Shapes[i].comp = icomp;
1136
1137 /*
1138 printf("i = %d, shape = %d, sign = %d, comp = %d\n", i, aRgn->Shapes[i].shape, aRgn->Shapes[i].sign, aRgn->Shapes[i].comp);
1139 */
1140
1141 }
1142
1143 return;
1144
1145 }
1146
1147 /*---------------------------------------------------------------------------*/
fits_setup_shape(RgnShape * newShape)1148 void fits_setup_shape ( RgnShape *newShape)
1149 {
1150 /* Perform some useful calculations now to speed up filter later */
1151
1152 double X, Y, R;
1153 double *coords;
1154 int i;
1155
1156 if ( newShape->shape == poly_rgn ) {
1157 coords = newShape->param.poly.Pts;
1158 } else {
1159 coords = newShape->param.gen.p;
1160 }
1161
1162 switch( newShape->shape ) {
1163 case circle_rgn:
1164 newShape->param.gen.a = coords[2] * coords[2];
1165 break;
1166 case annulus_rgn:
1167 newShape->param.gen.a = coords[2] * coords[2];
1168 newShape->param.gen.b = coords[3] * coords[3];
1169 break;
1170 case sector_rgn:
1171 while( coords[2]> 180.0 ) coords[2] -= 360.0;
1172 while( coords[2]<=-180.0 ) coords[2] += 360.0;
1173 while( coords[3]> 180.0 ) coords[3] -= 360.0;
1174 while( coords[3]<=-180.0 ) coords[3] += 360.0;
1175 break;
1176 case ellipse_rgn:
1177 newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
1178 newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
1179 break;
1180 case elliptannulus_rgn:
1181 newShape->param.gen.a = sin( myPI * (coords[6] / 180.0) );
1182 newShape->param.gen.b = cos( myPI * (coords[6] / 180.0) );
1183 newShape->param.gen.sinT = sin( myPI * (coords[7] / 180.0) );
1184 newShape->param.gen.cosT = cos( myPI * (coords[7] / 180.0) );
1185 break;
1186 case box_rgn:
1187 newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
1188 newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
1189 break;
1190 case boxannulus_rgn:
1191 newShape->param.gen.a = sin( myPI * (coords[6] / 180.0) );
1192 newShape->param.gen.b = cos( myPI * (coords[6] / 180.0) );
1193 newShape->param.gen.sinT = sin( myPI * (coords[7] / 180.0) );
1194 newShape->param.gen.cosT = cos( myPI * (coords[7] / 180.0) );
1195 break;
1196 case rectangle_rgn:
1197 newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
1198 newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
1199 X = 0.5 * ( coords[2]-coords[0] );
1200 Y = 0.5 * ( coords[3]-coords[1] );
1201 newShape->param.gen.a = fabs( X * newShape->param.gen.cosT
1202 + Y * newShape->param.gen.sinT );
1203 newShape->param.gen.b = fabs( Y * newShape->param.gen.cosT
1204 - X * newShape->param.gen.sinT );
1205 newShape->param.gen.p[5] = 0.5 * ( coords[2]+coords[0] );
1206 newShape->param.gen.p[6] = 0.5 * ( coords[3]+coords[1] );
1207 break;
1208 case diamond_rgn:
1209 newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
1210 newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
1211 break;
1212 case line_rgn:
1213 X = coords[2] - coords[0];
1214 Y = coords[3] - coords[1];
1215 R = sqrt( X*X + Y*Y );
1216 newShape->param.gen.sinT = ( R ? Y/R : 0.0 );
1217 newShape->param.gen.cosT = ( R ? X/R : 1.0 );
1218 newShape->param.gen.a = R + 0.5;
1219 break;
1220 case panda_rgn:
1221 while( coords[2]> 180.0 ) coords[2] -= 360.0;
1222 while( coords[2]<=-180.0 ) coords[2] += 360.0;
1223 while( coords[3]> 180.0 ) coords[3] -= 360.0;
1224 while( coords[3]<=-180.0 ) coords[3] += 360.0;
1225 newShape->param.gen.a = newShape->param.gen.p[5]*newShape->param.gen.p[5];
1226 newShape->param.gen.b = newShape->param.gen.p[6]*newShape->param.gen.p[6];
1227 break;
1228 case epanda_rgn:
1229 case bpanda_rgn:
1230 while( coords[2]> 180.0 ) coords[2] -= 360.0;
1231 while( coords[2]<=-180.0 ) coords[2] += 360.0;
1232 while( coords[3]> 180.0 ) coords[3] -= 360.0;
1233 while( coords[3]<=-180.0 ) coords[3] += 360.0;
1234 newShape->param.gen.sinT = sin( myPI * (coords[10] / 180.0) );
1235 newShape->param.gen.cosT = cos( myPI * (coords[10] / 180.0) );
1236 break;
1237 default:
1238 break;
1239 }
1240
1241 /* Set the xmin, xmax, ymin, ymax elements of the RgnShape structure */
1242
1243 /* For everything which has first two parameters as center position just */
1244 /* find a circle that encompasses the region and use it to set the */
1245 /* bounding box */
1246
1247 R = -1.0;
1248
1249 switch ( newShape->shape ) {
1250
1251 case circle_rgn:
1252 R = coords[2];
1253 break;
1254
1255 case annulus_rgn:
1256 R = coords[3];
1257 break;
1258
1259 case ellipse_rgn:
1260 if ( coords[2] > coords[3] ) {
1261 R = coords[2];
1262 } else {
1263 R = coords[3];
1264 }
1265 break;
1266
1267 case elliptannulus_rgn:
1268 if ( coords[4] > coords[5] ) {
1269 R = coords[4];
1270 } else {
1271 R = coords[5];
1272 }
1273 break;
1274
1275 case box_rgn:
1276 R = sqrt(coords[2]*coords[2]+
1277 coords[3]*coords[3])/2.0;
1278 break;
1279
1280 case boxannulus_rgn:
1281 R = sqrt(coords[4]*coords[5]+
1282 coords[4]*coords[5])/2.0;
1283 break;
1284
1285 case diamond_rgn:
1286 if ( coords[2] > coords[3] ) {
1287 R = coords[2]/2.0;
1288 } else {
1289 R = coords[3]/2.0;
1290 }
1291 break;
1292
1293 case point_rgn:
1294 R = 1.0;
1295 break;
1296
1297 case panda_rgn:
1298 R = coords[6];
1299 break;
1300
1301 case epanda_rgn:
1302 if ( coords[7] > coords[8] ) {
1303 R = coords[7];
1304 } else {
1305 R = coords[8];
1306 }
1307 break;
1308
1309 case bpanda_rgn:
1310 R = sqrt(coords[7]*coords[8]+
1311 coords[7]*coords[8])/2.0;
1312 break;
1313
1314 default:
1315 break;
1316 }
1317
1318 if ( R > 0.0 ) {
1319
1320 newShape->xmin = coords[0] - R;
1321 newShape->xmax = coords[0] + R;
1322 newShape->ymin = coords[1] - R;
1323 newShape->ymax = coords[1] + R;
1324
1325 return;
1326
1327 }
1328
1329 /* Now do the rest of the shapes that require individual methods */
1330
1331 switch ( newShape->shape ) {
1332
1333 case rectangle_rgn:
1334 R = sqrt((coords[5]-coords[0])*(coords[5]-coords[0])+
1335 (coords[6]-coords[1])*(coords[6]-coords[1]));
1336 newShape->xmin = coords[5] - R;
1337 newShape->xmax = coords[5] + R;
1338 newShape->ymin = coords[6] - R;
1339 newShape->ymax = coords[6] + R;
1340 break;
1341
1342 case poly_rgn:
1343 newShape->xmin = coords[0];
1344 newShape->xmax = coords[0];
1345 newShape->ymin = coords[1];
1346 newShape->ymax = coords[1];
1347 for( i=2; i < newShape->param.poly.nPts; ) {
1348 if( newShape->xmin > coords[i] ) /* Min X */
1349 newShape->xmin = coords[i];
1350 if( newShape->xmax < coords[i] ) /* Max X */
1351 newShape->xmax = coords[i];
1352 i++;
1353 if( newShape->ymin > coords[i] ) /* Min Y */
1354 newShape->ymin = coords[i];
1355 if( newShape->ymax < coords[i] ) /* Max Y */
1356 newShape->ymax = coords[i];
1357 i++;
1358 }
1359 break;
1360
1361 case line_rgn:
1362 if ( coords[0] > coords[2] ) {
1363 newShape->xmin = coords[2];
1364 newShape->xmax = coords[0];
1365 } else {
1366 newShape->xmin = coords[0];
1367 newShape->xmax = coords[2];
1368 }
1369 if ( coords[1] > coords[3] ) {
1370 newShape->ymin = coords[3];
1371 newShape->ymax = coords[1];
1372 } else {
1373 newShape->ymin = coords[1];
1374 newShape->ymax = coords[3];
1375 }
1376
1377 break;
1378
1379 /* sector doesn't have min and max so indicate by setting max < min */
1380
1381 case sector_rgn:
1382 newShape->xmin = 1.0;
1383 newShape->xmax = -1.0;
1384 newShape->ymin = 1.0;
1385 newShape->ymax = -1.0;
1386 break;
1387
1388 default:
1389 break;
1390 }
1391
1392 return;
1393
1394 }
1395
1396 /*---------------------------------------------------------------------------*/
fits_read_fits_region(fitsfile * fptr,WCSdata * wcs,SAORegion ** Rgn,int * status)1397 int fits_read_fits_region ( fitsfile *fptr,
1398 WCSdata *wcs,
1399 SAORegion **Rgn,
1400 int *status)
1401 /* Read regions from a FITS region extension and return the information */
1402 /* in the "SAORegion" structure. If it is nonNULL, use wcs to convert the */
1403 /* region coordinates to pixels. Return an error if region is in degrees */
1404 /* but no WCS data is provided. */
1405 /*---------------------------------------------------------------------------*/
1406 {
1407
1408 int i, j, icol[6], idum, anynul, npos;
1409 int dotransform, got_component = 1, tstatus;
1410 long icsize[6];
1411 double X, Y, Theta, Xsave = 0, Ysave = 0, Xpos, Ypos;
1412 double *coords;
1413 char *cvalue, *cvalue2;
1414 char comment[FLEN_COMMENT];
1415 char colname[6][FLEN_VALUE] = {"X", "Y", "SHAPE", "R", "ROTANG", "COMPONENT"};
1416 char shapename[17][FLEN_VALUE] = {"POINT","CIRCLE","ELLIPSE","ANNULUS",
1417 "ELLIPTANNULUS","BOX","ROTBOX","BOXANNULUS",
1418 "RECTANGLE","ROTRECTANGLE","POLYGON","PIE",
1419 "SECTOR","DIAMOND","RHOMBUS","ROTDIAMOND",
1420 "ROTRHOMBUS"};
1421 int shapetype[17] = {point_rgn, circle_rgn, ellipse_rgn, annulus_rgn,
1422 elliptannulus_rgn, box_rgn, box_rgn, boxannulus_rgn,
1423 rectangle_rgn, rectangle_rgn, poly_rgn, sector_rgn,
1424 sector_rgn, diamond_rgn, diamond_rgn, diamond_rgn,
1425 diamond_rgn};
1426 SAORegion *aRgn;
1427 RgnShape *newShape;
1428 WCSdata *regwcs = 0;
1429
1430 if ( *status ) return( *status );
1431
1432 aRgn = (SAORegion *)malloc( sizeof(SAORegion) );
1433 if( ! aRgn ) {
1434 ffpmsg("Couldn't allocate memory to hold Region file contents.");
1435 return(*status = MEMORY_ALLOCATION );
1436 }
1437 aRgn->nShapes = 0;
1438 aRgn->Shapes = NULL;
1439 if( wcs && wcs->exists )
1440 aRgn->wcs = *wcs;
1441 else
1442 aRgn->wcs.exists = 0;
1443
1444 /* See if we are already positioned to a region extension, else */
1445 /* move to the REGION extension (file is already open). */
1446
1447 tstatus = 0;
1448 for (i=0; i<5; i++) {
1449 ffgcno(fptr, CASEINSEN, colname[i], &icol[i], &tstatus);
1450 }
1451
1452 if (tstatus) {
1453 /* couldn't find the required columns, so search for "REGION" extension */
1454 if ( ffmnhd(fptr, BINARY_TBL, "REGION", 1, status) ) {
1455 ffpmsg("Could not move to REGION extension.");
1456 goto error;
1457 }
1458 }
1459
1460 /* get the number of shapes and allocate memory */
1461
1462 if ( ffgky(fptr, TINT, "NAXIS2", &aRgn->nShapes, comment, status) ) {
1463 ffpmsg("Could not read NAXIS2 keyword.");
1464 goto error;
1465 }
1466
1467 aRgn->Shapes = (RgnShape *) malloc(aRgn->nShapes * sizeof(RgnShape));
1468 if ( !aRgn->Shapes ) {
1469 ffpmsg( "Failed to allocate memory for Region data");
1470 *status = MEMORY_ALLOCATION;
1471 goto error;
1472 }
1473
1474 /* get the required column numbers */
1475
1476 for (i=0; i<5; i++) {
1477 if ( ffgcno(fptr, CASEINSEN, colname[i], &icol[i], status) ) {
1478 ffpmsg("Could not find column.");
1479 goto error;
1480 }
1481 }
1482
1483 /* try to get the optional column numbers */
1484
1485 if ( ffgcno(fptr, CASEINSEN, colname[5], &icol[5], status) ) {
1486 got_component = 0;
1487 }
1488
1489 /* if there was input WCS then read the WCS info for the region in case they */
1490 /* are different and we have to transform */
1491
1492 dotransform = 0;
1493 if ( aRgn->wcs.exists ) {
1494 regwcs = (WCSdata *) malloc ( sizeof(WCSdata) );
1495 if ( !regwcs ) {
1496 ffpmsg( "Failed to allocate memory for Region WCS data");
1497 *status = MEMORY_ALLOCATION;
1498 goto error;
1499 }
1500
1501 regwcs->exists = 1;
1502 if ( ffgtcs(fptr, icol[0], icol[1], ®wcs->xrefval, ®wcs->yrefval,
1503 ®wcs->xrefpix, ®wcs->yrefpix, ®wcs->xinc, ®wcs->yinc,
1504 ®wcs->rot, regwcs->type, status) ) {
1505 regwcs->exists = 0;
1506 *status = 0;
1507 }
1508
1509 if ( regwcs->exists && wcs->exists ) {
1510 if ( fabs(regwcs->xrefval-wcs->xrefval) > 1.0e-6 ||
1511 fabs(regwcs->yrefval-wcs->yrefval) > 1.0e-6 ||
1512 fabs(regwcs->xrefpix-wcs->xrefpix) > 1.0e-6 ||
1513 fabs(regwcs->yrefpix-wcs->yrefpix) > 1.0e-6 ||
1514 fabs(regwcs->xinc-wcs->xinc) > 1.0e-6 ||
1515 fabs(regwcs->yinc-wcs->yinc) > 1.0e-6 ||
1516 fabs(regwcs->rot-wcs->rot) > 1.0e-6 ||
1517 !strcmp(regwcs->type,wcs->type) ) dotransform = 1;
1518 }
1519 }
1520
1521 /* get the sizes of the X, Y, R, and ROTANG vectors */
1522
1523 for (i=0; i<6; i++) {
1524 if ( ffgtdm(fptr, icol[i], 1, &idum, &icsize[i], status) ) {
1525 ffpmsg("Could not find vector size of column.");
1526 goto error;
1527 }
1528 }
1529
1530 cvalue = (char *) malloc ((FLEN_VALUE+1)*sizeof(char));
1531
1532 /* loop over the shapes - note 1-based counting for rows in FITS files */
1533
1534 for (i=1; i<=aRgn->nShapes; i++) {
1535
1536 newShape = &aRgn->Shapes[i-1];
1537 for (j=0; j<8; j++) newShape->param.gen.p[j] = 0.0;
1538 newShape->param.gen.a = 0.0;
1539 newShape->param.gen.b = 0.0;
1540 newShape->param.gen.sinT = 0.0;
1541 newShape->param.gen.cosT = 0.0;
1542
1543 /* get the shape */
1544
1545 if ( ffgcvs(fptr, icol[2], i, 1, 1, " ", &cvalue, &anynul, status) ) {
1546 ffpmsg("Could not read shape.");
1547 goto error;
1548 }
1549
1550 /* set include or exclude */
1551
1552 newShape->sign = 1;
1553 cvalue2 = cvalue;
1554 if ( !strncmp(cvalue,"!",1) ) {
1555 newShape->sign = 0;
1556 cvalue2++;
1557 }
1558
1559 /* set the shape type */
1560
1561 for (j=0; j<17; j++) {
1562 if ( !strcmp(cvalue2, shapename[j]) ) newShape->shape = shapetype[j];
1563 }
1564
1565 /* allocate memory for polygon case and set coords pointer */
1566
1567 if ( newShape->shape == poly_rgn ) {
1568 newShape->param.poly.Pts = (double *) calloc (2*icsize[0], sizeof(double));
1569 if ( !newShape->param.poly.Pts ) {
1570 ffpmsg("Could not allocate memory to hold polygon parameters" );
1571 *status = MEMORY_ALLOCATION;
1572 goto error;
1573 }
1574 newShape->param.poly.nPts = 2*icsize[0];
1575 coords = newShape->param.poly.Pts;
1576 } else {
1577 coords = newShape->param.gen.p;
1578 }
1579
1580
1581 /* read X and Y. Polygon and Rectangle require special cases */
1582
1583 npos = 1;
1584 if ( newShape->shape == poly_rgn ) npos = newShape->param.poly.nPts/2;
1585 if ( newShape->shape == rectangle_rgn ) npos = 2;
1586
1587 for (j=0; j<npos; j++) {
1588 if ( ffgcvd(fptr, icol[0], i, j+1, 1, DOUBLENULLVALUE, coords, &anynul, status) ) {
1589 ffpmsg("Failed to read X column for polygon region");
1590 goto error;
1591 }
1592 if (*coords == DOUBLENULLVALUE) { /* check for null value end of array marker */
1593 npos = j;
1594 newShape->param.poly.nPts = npos * 2;
1595 break;
1596 }
1597 coords++;
1598
1599 if ( ffgcvd(fptr, icol[1], i, j+1, 1, DOUBLENULLVALUE, coords, &anynul, status) ) {
1600 ffpmsg("Failed to read Y column for polygon region");
1601 goto error;
1602 }
1603 if (*coords == DOUBLENULLVALUE) { /* check for null value end of array marker */
1604 npos = j;
1605 newShape->param.poly.nPts = npos * 2;
1606 coords--;
1607 break;
1608 }
1609 coords++;
1610
1611 if (j == 0) { /* save the first X and Y coordinate */
1612 Xsave = *(coords - 2);
1613 Ysave = *(coords - 1);
1614 } else if ((Xsave == *(coords - 2)) && (Ysave == *(coords - 1)) ) {
1615 /* if point has same coordinate as first point, this marks the end of the array */
1616 npos = j + 1;
1617 newShape->param.poly.nPts = npos * 2;
1618 break;
1619 }
1620 }
1621
1622 /* transform positions if the region and input wcs differ */
1623
1624 if ( dotransform ) {
1625
1626 coords -= npos*2;
1627 Xsave = coords[0];
1628 Ysave = coords[1];
1629 for (j=0; j<npos; j++) {
1630 ffwldp(coords[2*j], coords[2*j+1], regwcs->xrefval, regwcs->yrefval, regwcs->xrefpix,
1631 regwcs->yrefpix, regwcs->xinc, regwcs->yinc, regwcs->rot,
1632 regwcs->type, &Xpos, &Ypos, status);
1633 ffxypx(Xpos, Ypos, wcs->xrefval, wcs->yrefval, wcs->xrefpix,
1634 wcs->yrefpix, wcs->xinc, wcs->yinc, wcs->rot,
1635 wcs->type, &coords[2*j], &coords[2*j+1], status);
1636 if ( *status ) {
1637 ffpmsg("Failed to transform coordinates");
1638 goto error;
1639 }
1640 }
1641 coords += npos*2;
1642 }
1643
1644 /* read R. Circle requires one number; Box, Diamond, Ellipse, Annulus, Sector
1645 and Panda two; Boxannulus and Elliptannulus four; Point, Rectangle and
1646 Polygon none. */
1647
1648 npos = 0;
1649 switch ( newShape->shape ) {
1650 case circle_rgn:
1651 npos = 1;
1652 break;
1653 case box_rgn:
1654 case diamond_rgn:
1655 case ellipse_rgn:
1656 case annulus_rgn:
1657 case sector_rgn:
1658 npos = 2;
1659 break;
1660 case boxannulus_rgn:
1661 case elliptannulus_rgn:
1662 npos = 4;
1663 break;
1664 default:
1665 break;
1666 }
1667
1668 if ( npos > 0 ) {
1669 if ( ffgcvd(fptr, icol[3], i, 1, npos, 0.0, coords, &anynul, status) ) {
1670 ffpmsg("Failed to read R column for region");
1671 goto error;
1672 }
1673
1674 /* transform lengths if the region and input wcs differ */
1675
1676 if ( dotransform ) {
1677 for (j=0; j<npos; j++) {
1678 Y = Ysave + (*coords);
1679 X = Xsave;
1680 ffwldp(X, Y, regwcs->xrefval, regwcs->yrefval, regwcs->xrefpix,
1681 regwcs->yrefpix, regwcs->xinc, regwcs->yinc, regwcs->rot,
1682 regwcs->type, &Xpos, &Ypos, status);
1683 ffxypx(Xpos, Ypos, wcs->xrefval, wcs->yrefval, wcs->xrefpix,
1684 wcs->yrefpix, wcs->xinc, wcs->yinc, wcs->rot,
1685 wcs->type, &X, &Y, status);
1686 if ( *status ) {
1687 ffpmsg("Failed to transform coordinates");
1688 goto error;
1689 }
1690 *(coords++) = sqrt(pow(X-newShape->param.gen.p[0],2)+pow(Y-newShape->param.gen.p[1],2));
1691 }
1692 } else {
1693 coords += npos;
1694 }
1695 }
1696
1697 /* read ROTANG. Requires two values for Boxannulus, Elliptannulus, Sector,
1698 Panda; one for Box, Diamond, Ellipse; and none for Circle, Point, Annulus,
1699 Rectangle, Polygon */
1700
1701 npos = 0;
1702 switch ( newShape->shape ) {
1703 case box_rgn:
1704 case diamond_rgn:
1705 case ellipse_rgn:
1706 npos = 1;
1707 break;
1708 case boxannulus_rgn:
1709 case elliptannulus_rgn:
1710 case sector_rgn:
1711 npos = 2;
1712 break;
1713 default:
1714 break;
1715 }
1716
1717 if ( npos > 0 ) {
1718 if ( ffgcvd(fptr, icol[4], i, 1, npos, 0.0, coords, &anynul, status) ) {
1719 ffpmsg("Failed to read ROTANG column for region");
1720 goto error;
1721 }
1722
1723 /* transform angles if the region and input wcs differ */
1724
1725 if ( dotransform ) {
1726 Theta = (wcs->rot) - (regwcs->rot);
1727 for (j=0; j<npos; j++) *(coords++) += Theta;
1728 } else {
1729 coords += npos;
1730 }
1731 }
1732
1733 /* read the component number */
1734
1735 if (got_component) {
1736 if ( ffgcv(fptr, TINT, icol[5], i, 1, 1, 0, &newShape->comp, &anynul, status) ) {
1737 ffpmsg("Failed to read COMPONENT column for region");
1738 goto error;
1739 }
1740 } else {
1741 newShape->comp = 1;
1742 }
1743
1744
1745 /* do some precalculations to speed up tests */
1746
1747 fits_setup_shape(newShape);
1748
1749 /* end loop over shapes */
1750
1751 }
1752
1753 error:
1754
1755 if( *status )
1756 fits_free_region( aRgn );
1757 else
1758 *Rgn = aRgn;
1759
1760 ffclos(fptr, status);
1761
1762 return( *status );
1763 }
1764
1765