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