« back to all changes in this revision
Viewing changes to region.c
- browse files at revision 1
- compare with another revision
- download diff
- download tarball
- view history from revision 1
- Committer: Bazaar Package Importer
- Author(s): Gopal Narayanan
- Date: 2002-02-26 11:27:29 UTC
- Revision ID: james.westby@ubuntu.com-20020226112729-3q2o993rhh81ipp4
Tags: upstream-2.401
Import upstream version 2.401
- files added:
added
removed
region.c
1
#include <stdio.h>
2
#include <stdlib.h>
3
#include <string.h>
4
#include <math.h>
5
#include "fitsio2.h"
6
#include "region.h"
7
8
static int Pt_in_Poly( double x, double y, int nPts, double *Pts );
9
10
/*---------------------------------------------------------------------------*/
11
int ffrrgn( const char *filename,
12
WCSdata *wcs,
13
SAORegion **Rgn,
14
int *status )
15
/* Read regions from a SAO-style 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
char *currLine;
22
char *namePtr, *paramPtr, *currLoc;
23
char *pX, *pY, *endp;
24
long allocLen, lineLen, hh, mm, dd;
25
double *coords = 0, X, Y, R, x, y, ss, xsave= 0., ysave= 0.;
26
int nParams, nCoords, negdec;
27
int i, done;
28
FILE *rgnFile;
29
coordFmt cFmt;
30
SAORegion *aRgn;
31
RgnShape *newShape, *tmpShape;
32
33
if( *status ) return( *status );
34
35
aRgn = (SAORegion *)malloc( sizeof(SAORegion) );
36
if( ! aRgn ) {
37
ffpmsg("Couldn't allocate memory to hold Region file contents.");
38
return(*status = MEMORY_ALLOCATION );
39
}
40
aRgn->nShapes = 0;
41
aRgn->Shapes = NULL;
42
if( wcs && wcs->exists )
43
aRgn->wcs = *wcs;
44
else
45
aRgn->wcs.exists = 0;
46
47
cFmt = pixel_fmt; /* set default format */
48
49
/* Allocate Line Buffer */
50
51
allocLen = 512;
52
currLine = (char *)malloc( allocLen * sizeof(char) );
53
if( !currLine ) {
54
free( aRgn );
55
ffpmsg("Couldn't allocate memory to hold Region file contents.");
56
return(*status = MEMORY_ALLOCATION );
57
}
58
59
/* Open Region File */
60
61
if( (rgnFile = fopen( filename, "r" ))==NULL ) {
62
sprintf(currLine,"Could not open Region file %s.",filename);
63
ffpmsg( currLine );
64
free( currLine );
65
free( aRgn );
66
return( *status = FILE_NOT_OPENED );
67
}
68
69
/* Read in file, line by line */
70
71
while( fgets(currLine,allocLen,rgnFile) != NULL ) {
72
73
/* Make sure we have a full line of text */
74
75
lineLen = strlen(currLine);
76
while( lineLen==allocLen-1 && currLine[lineLen-1]!='\n' ) {
77
currLoc = (char *)realloc( currLine, 2 * allocLen * sizeof(char) );
78
if( !currLoc ) {
79
ffpmsg("Couldn't allocate memory to hold Region file contents.");
80
*status = MEMORY_ALLOCATION;
81
goto error;
82
} else {
83
currLine = currLoc;
84
}
85
fgets( currLine+lineLen, allocLen+1, rgnFile );
86
allocLen += allocLen;
87
lineLen += strlen(currLine+lineLen);
88
}
89
90
currLoc = currLine;
91
if( *currLoc == '#' ) {
92
93
/* Look to see if it is followed by a format statement... */
94
/* if not skip line */
95
96
currLoc++;
97
while( *currLoc==' ' ) currLoc++;
98
if( !strncasecmp( currLoc, "format:", 7 ) ) {
99
if( aRgn->nShapes ) {
100
ffpmsg("Format code encountered after reading 1 or more shapes.");
101
*status = PARSE_SYNTAX_ERR;
102
goto error;
103
}
104
currLoc += 7;
105
while( *currLoc==' ' ) currLoc++;
106
if( !strncasecmp( currLoc, "pixel", 5 ) ) {
107
cFmt = pixel_fmt;
108
} else if( !strncasecmp( currLoc, "degree", 6 ) ) {
109
cFmt = degree_fmt;
110
} else if( !strncasecmp( currLoc, "hhmmss", 6 ) ) {
111
cFmt = hhmmss_fmt;
112
} else if( !strncasecmp( currLoc, "hms", 3 ) ) {
113
cFmt = hhmmss_fmt;
114
} else {
115
ffpmsg("Unknown format code encountered in region file.");
116
*status = PARSE_SYNTAX_ERR;
117
goto error;
118
}
119
}
120
121
} else if( !strncasecmp( currLoc, "glob", 4 ) ) {
122
/* skip lines that begin with the word 'global' */
123
124
} else {
125
126
while( *currLoc != '\0' ) {
127
128
namePtr = currLoc;
129
paramPtr = NULL;
130
nParams = 1;
131
132
/* Search for closing parenthesis */
133
134
done = 0;
135
while( !done && !*status && *currLoc ) {
136
switch (*currLoc) {
137
case '(':
138
*currLoc = '\0';
139
currLoc++;
140
if( paramPtr ) /* Can't have two '(' in a region! */
141
*status = 1;
142
else
143
paramPtr = currLoc;
144
break;
145
case ')':
146
*currLoc = '\0';
147
currLoc++;
148
if( !paramPtr ) /* Can't have a ')' without a '(' first */
149
*status = 1;
150
else
151
done = 1;
152
break;
153
case '#':
154
case '\n':
155
*currLoc = '\0';
156
if( !paramPtr ) /* Allow for a blank line */
157
done = 1;
158
break;
159
case ':':
160
currLoc++;
161
cFmt = hhmmss_fmt;
162
break;
163
case 'd':
164
currLoc++;
165
cFmt = degree_fmt;
166
break;
167
case ',':
168
nParams++; /* Fall through to default */
169
default:
170
currLoc++;
171
break;
172
}
173
}
174
if( *status || !done ) {
175
ffpmsg( "Error reading Region file" );
176
*status = PARSE_SYNTAX_ERR;
177
goto error;
178
}
179
180
/* Skip white space in region name */
181
182
while( *namePtr==' ' ) namePtr++;
183
184
/* Was this a blank line? Or the end of the current one */
185
186
if( ! *namePtr && ! paramPtr ) continue;
187
188
/**************************************************/
189
/* We've apparently found a region... Set it up */
190
/**************************************************/
191
192
if( !(aRgn->nShapes % 10) ) {
193
if( aRgn->Shapes )
194
tmpShape = (RgnShape *)realloc( aRgn->Shapes,
195
(10+aRgn->nShapes)
196
* sizeof(RgnShape) );
197
else
198
tmpShape = (RgnShape *) malloc( 10 * sizeof(RgnShape) );
199
if( tmpShape ) {
200
aRgn->Shapes = tmpShape;
201
} else {
202
ffpmsg( "Failed to allocate memory for Region data");
203
*status = MEMORY_ALLOCATION;
204
goto error;
205
}
206
207
}
208
newShape = &aRgn->Shapes[aRgn->nShapes++];
209
newShape->sign = 1;
210
newShape->shape = point_rgn;
211
212
/* Check for format code at beginning of the line */
213
214
if( !strncasecmp( namePtr, "image;", 6 ) ) {
215
namePtr += 6;
216
cFmt = pixel_fmt;
217
} else if( !strncasecmp( namePtr, "fk4;", 4 ) ) {
218
namePtr += 4;
219
cFmt = degree_fmt;
220
} else if( !strncasecmp( namePtr, "fk5;", 4 ) ) {
221
namePtr += 4;
222
cFmt = degree_fmt;
223
} else if( !strncasecmp( namePtr, "icrs;", 5 ) ) {
224
namePtr += 5;
225
cFmt = degree_fmt;
226
} else if( !strncasecmp( namePtr, "fk5", 3 ) ) {
227
cFmt = degree_fmt;
228
continue; /* supports POW region file format */
229
} else if( !strncasecmp( namePtr, "fk4", 3 ) ) {
230
cFmt = degree_fmt;
231
continue; /* supports POW region file format */
232
} else if( !strncasecmp( namePtr, "icrs", 4 ) ) {
233
cFmt = degree_fmt;
234
continue; /* supports POW region file format */
235
} else if( !strncasecmp( namePtr, "image", 5 ) ) {
236
cFmt = pixel_fmt;
237
continue; /* supports POW region file format */
238
} else if( !strncasecmp( namePtr, "galactic;", 9 ) ) {
239
ffpmsg( "Galactic region coordinates not supported" );
240
ffpmsg( namePtr );
241
*status = PARSE_SYNTAX_ERR;
242
goto error;
243
} else if( !strncasecmp( namePtr, "ecliptic;", 9 ) ) {
244
ffpmsg( "ecliptic region coordinates not supported" );
245
ffpmsg( namePtr );
246
*status = PARSE_SYNTAX_ERR;
247
goto error;
248
}
249
250
while( *namePtr==' ' ) namePtr++;
251
252
/* Check for the shape's sign */
253
254
if( *namePtr=='+' ) {
255
namePtr++;
256
} else if( *namePtr=='-' ) {
257
namePtr++;
258
newShape->sign = 0;
259
}
260
261
/* Skip white space in region name */
262
263
while( *namePtr==' ' ) namePtr++;
264
if( *namePtr=='\0' ) {
265
ffpmsg( "Error reading Region file" );
266
*status = PARSE_SYNTAX_ERR;
267
goto error;
268
}
269
lineLen = strlen( namePtr ) - 1;
270
while( namePtr[lineLen]==' ' ) namePtr[lineLen--] = '\0';
271
272
/* Now identify the region */
273
274
if( !strcasecmp( namePtr, "circle" ) ) {
275
newShape->shape = circle_rgn;
276
if( nParams != 3 )
277
*status = PARSE_SYNTAX_ERR;
278
nCoords = 2;
279
} else if( !strcasecmp( namePtr, "annulus" ) ) {
280
newShape->shape = annulus_rgn;
281
if( nParams != 4 )
282
*status = PARSE_SYNTAX_ERR;
283
nCoords = 2;
284
} else if( !strcasecmp( namePtr, "ellipse" ) ) {
285
newShape->shape = ellipse_rgn;
286
if( nParams < 4 || nParams > 5 )
287
*status = PARSE_SYNTAX_ERR;
288
newShape->param.gen.p[4] = 0.0;
289
nCoords = 2;
290
} else if( !strcasecmp( namePtr, "elliptannulus" ) ) {
291
newShape->shape = elliptannulus_rgn;
292
if( !( nParams==8 || nParams==6 ) )
293
*status = PARSE_SYNTAX_ERR;
294
newShape->param.gen.p[6] = 0.0;
295
newShape->param.gen.p[7] = 0.0;
296
nCoords = 2;
297
} else if( !strcasecmp( namePtr, "box" )
298
|| !strcasecmp( namePtr, "rotbox" ) ) {
299
newShape->shape = box_rgn;
300
if( nParams < 4 || nParams > 5 )
301
*status = PARSE_SYNTAX_ERR;
302
newShape->param.gen.p[4] = 0.0;
303
nCoords = 2;
304
} else if( !strcasecmp( namePtr, "rectangle" )
305
|| !strcasecmp( namePtr, "rotrectangle" ) ) {
306
newShape->shape = rectangle_rgn;
307
if( nParams < 4 || nParams > 5 )
308
*status = PARSE_SYNTAX_ERR;
309
newShape->param.gen.p[4] = 0.0;
310
nCoords = 4;
311
} else if( !strcasecmp( namePtr, "diamond" )
312
|| !strcasecmp( namePtr, "rotdiamond" )
313
|| !strcasecmp( namePtr, "rhombus" )
314
|| !strcasecmp( namePtr, "rotrhombus" ) ) {
315
newShape->shape = diamond_rgn;
316
if( nParams < 4 || nParams > 5 )
317
*status = PARSE_SYNTAX_ERR;
318
newShape->param.gen.p[4] = 0.0;
319
nCoords = 2;
320
} else if( !strcasecmp( namePtr, "sector" )
321
|| !strcasecmp( namePtr, "pie" ) ) {
322
newShape->shape = sector_rgn;
323
if( nParams != 4 )
324
*status = PARSE_SYNTAX_ERR;
325
nCoords = 2;
326
} else if( !strcasecmp( namePtr, "point" ) ) {
327
newShape->shape = point_rgn;
328
if( nParams != 2 )
329
*status = PARSE_SYNTAX_ERR;
330
nCoords = 2;
331
} else if( !strcasecmp( namePtr, "line" ) ) {
332
newShape->shape = line_rgn;
333
if( nParams != 4 )
334
*status = PARSE_SYNTAX_ERR;
335
nCoords = 4;
336
} else if( !strcasecmp( namePtr, "polygon" ) ) {
337
newShape->shape = poly_rgn;
338
if( nParams < 6 || (nParams&1) )
339
*status = PARSE_SYNTAX_ERR;
340
nCoords = nParams;
341
} else {
342
ffpmsg( "Unrecognized region found in region file:" );
343
ffpmsg( namePtr );
344
*status = PARSE_SYNTAX_ERR;
345
goto error;
346
}
347
if( *status ) {
348
ffpmsg( "Wrong number of parameters found for region" );
349
ffpmsg( namePtr );
350
goto error;
351
}
352
353
/* Parse Parameter string... convert to pixels if necessary */
354
355
if( newShape->shape==poly_rgn ) {
356
newShape->param.poly.Pts = (double *)malloc( nParams
357
* sizeof(double) );
358
if( !newShape->param.poly.Pts ) {
359
ffpmsg(
360
"Could not allocate memory to hold polygon parameters" );
361
*status = MEMORY_ALLOCATION;
362
goto error;
363
}
364
newShape->param.poly.nPts = nParams;
365
coords = newShape->param.poly.Pts;
366
} else
367
coords = newShape->param.gen.p;
368
369
/* Parse the initial "WCS?" coordinates */
370
for( i=0; i<nCoords; i+=2 ) {
371
pX = paramPtr;
372
while( *paramPtr!=',' ) paramPtr++;
373
*(paramPtr++) = '\0';
374
375
pY = paramPtr;
376
while( *paramPtr!=',' && *paramPtr != '\0' ) paramPtr++;
377
*(paramPtr++) = '\0';
378
379
if( strchr(pX, ':' ) ) {
380
/* Read in special format & convert to decimal degrees */
381
cFmt = hhmmss_fmt;
382
mm = 0;
383
ss = 0.;
384
hh = strtol(pX, &endp, 10);
385
if (endp && *endp==':') {
386
pX = endp + 1;
387
mm = strtol(pX, &endp, 10);
388
if (endp && *endp==':') {
389
pX = endp + 1;
390
ss = atof( pX );
391
}
392
}
393
X = 15. * (hh + mm/60. + ss/3600.); /* convert to degrees */
394
395
mm = 0;
396
ss = 0.;
397
negdec = 0;
398
399
while( *pY==' ' ) pY++;
400
if (*pY=='-') {
401
negdec = 1;
402
pY++;
403
}
404
dd = strtol(pY, &endp, 10);
405
if (endp && *endp==':') {
406
pY = endp + 1;
407
mm = strtol(pY, &endp, 10);
408
if (endp && *endp==':') {
409
pY = endp + 1;
410
ss = atof( pY );
411
}
412
}
413
if (negdec)
414
Y = -dd - mm/60. - ss/3600.; /* convert to degrees */
415
else
416
Y = dd + mm/60. + ss/3600.;
417
418
} else {
419
X = atof( pX );
420
Y = atof( pY );
421
}
422
if (i==0) { /* save 1st coord. in case needed later */
423
xsave = X;
424
ysave = Y;
425
}
426
427
if( cFmt!=pixel_fmt ) {
428
/* Convert to pixels */
429
if( wcs==NULL || ! wcs->exists ) {
430
ffpmsg("WCS information needed to convert region coordinates.");
431
*status = NO_WCS_KEY;
432
goto error;
433
}
434
435
if( ffxypx( X, Y, wcs->xrefval, wcs->yrefval,
436
wcs->xrefpix, wcs->yrefpix,
437
wcs->xinc, wcs->yinc,
438
wcs->rot, wcs->type,
439
&x, &y, status ) ) {
440
ffpmsg("Error converting region to pixel coordinates.");
441
goto error;
442
}
443
X = x; Y = y;
444
}
445
coords[i] = X;
446
coords[i+1] = Y;
447
}
448
449
/* Read in remaining parameters... */
450
451
for( ; i<nParams; i++ ) {
452
pX = paramPtr;
453
while( *paramPtr!=',' && *paramPtr != '\0' ) paramPtr++;
454
*(paramPtr++) = '\0';
455
coords[i] = strtod( pX, &endp );
456
457
if (endp && *endp=='"') {
458
/* parameter given in arcsec so convert to pixels. */
459
/* Increment first Y coordinate by this amount then calc */
460
/* the distance in pixels from the original coordinate. */
461
/* NOTE: This assumes the pixels are square!! */
462
if (ysave < 0.)
463
Y = ysave + coords[i]/3600.; /* don't exceed -90 */
464
else
465
Y = ysave - coords[i]/3600.; /* don't exceed +90 */
466
467
X = xsave;
468
if( ffxypx( X, Y, wcs->xrefval, wcs->yrefval,
469
wcs->xrefpix, wcs->yrefpix,
470
wcs->xinc, wcs->yinc,
471
wcs->rot, wcs->type,
472
&x, &y, status ) ) {
473
ffpmsg("Error converting region to pixel coordinates.");
474
goto error;
475
}
476
477
coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );
478
479
} else if (endp && *endp=='\'') {
480
/* parameter given in arcmin so convert to pixels. */
481
/* Increment first Y coordinate by this amount, then calc */
482
/* the distance in pixels from the original coordinate. */
483
/* NOTE: This assumes the pixels are square!! */
484
if (ysave < 0.)
485
Y = ysave + coords[i]/60.; /* don't exceed -90 */
486
else
487
Y = ysave - coords[i]/60.; /* don't exceed +90 */
488
489
X = xsave;
490
if( ffxypx( X, Y, wcs->xrefval, wcs->yrefval,
491
wcs->xrefpix, wcs->yrefpix,
492
wcs->xinc, wcs->yinc,
493
wcs->rot, wcs->type,
494
&x, &y, status ) ) {
495
ffpmsg("Error converting region to pixel coordinates.");
496
goto error;
497
}
498
499
coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );
500
501
} else if (endp && *endp=='d') {
502
/* parameter given in degrees so convert to pixels. */
503
/* Increment first Y coordinate by this amount, then calc */
504
/* the distance in pixels from the original coordinate. */
505
/* NOTE: This assumes the pixels are square!! */
506
if (ysave < 0.)
507
Y = ysave + coords[i]; /* don't exceed -90 */
508
else
509
Y = ysave - coords[i]; /* don't exceed +90 */
510
511
X = xsave;
512
if( ffxypx( X, Y, wcs->xrefval, wcs->yrefval,
513
wcs->xrefpix, wcs->yrefpix,
514
wcs->xinc, wcs->yinc,
515
wcs->rot, wcs->type,
516
&x, &y, status ) ) {
517
ffpmsg("Error converting region to pixel coordinates.");
518
goto error;
519
}
520
521
coords[i] = sqrt( pow(x-coords[0],2) + pow(y-coords[1],2) );
522
}
523
}
524
525
/* Perform some useful calculations now to speed up filter later */
526
527
switch( newShape->shape ) {
528
case circle_rgn:
529
newShape->param.gen.a = coords[2] * coords[2];
530
break;
531
case annulus_rgn:
532
newShape->param.gen.a = coords[2] * coords[2];
533
newShape->param.gen.b = coords[3] * coords[3];
534
break;
535
case sector_rgn:
536
while( coords[2]> 180.0 ) coords[2] -= 360.0;
537
while( coords[2]<=-180.0 ) coords[2] += 360.0;
538
while( coords[3]> 180.0 ) coords[3] -= 360.0;
539
while( coords[3]<=-180.0 ) coords[3] += 360.0;
540
break;
541
case ellipse_rgn:
542
newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
543
newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
544
break;
545
case elliptannulus_rgn:
546
newShape->param.gen.a = sin( myPI * (coords[6] / 180.0) );
547
newShape->param.gen.b = cos( myPI * (coords[6] / 180.0) );
548
newShape->param.gen.sinT = sin( myPI * (coords[7] / 180.0) );
549
newShape->param.gen.cosT = cos( myPI * (coords[7] / 180.0) );
550
break;
551
case box_rgn:
552
newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
553
newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
554
break;
555
case rectangle_rgn:
556
newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
557
newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
558
X = 0.5 * ( coords[2]-coords[0] );
559
Y = 0.5 * ( coords[3]-coords[1] );
560
newShape->param.gen.a = fabs( X * newShape->param.gen.cosT
561
+ Y * newShape->param.gen.sinT );
562
newShape->param.gen.b = fabs( Y * newShape->param.gen.cosT
563
- X * newShape->param.gen.sinT );
564
newShape->param.gen.p[5] = 0.5 * ( coords[2]+coords[0] );
565
newShape->param.gen.p[6] = 0.5 * ( coords[3]+coords[1] );
566
break;
567
case diamond_rgn:
568
newShape->param.gen.sinT = sin( myPI * (coords[4] / 180.0) );
569
newShape->param.gen.cosT = cos( myPI * (coords[4] / 180.0) );
570
break;
571
case line_rgn:
572
X = coords[2] - coords[0];
573
Y = coords[3] - coords[1];
574
R = sqrt( X*X + Y*Y );
575
newShape->param.gen.sinT = ( R ? Y/R : 0.0 );
576
newShape->param.gen.cosT = ( R ? X/R : 1.0 );
577
newShape->param.gen.a = R + 0.5;
578
break;
579
case point_rgn:
580
break;
581
case poly_rgn:
582
/* Find bounding box */
583
newShape->param.poly.xmin = coords[0];
584
newShape->param.poly.xmax = coords[0];
585
newShape->param.poly.ymin = coords[1];
586
newShape->param.poly.ymax = coords[1];
587
for( i=2; i<nParams; ) {
588
if( newShape->param.poly.xmin > coords[i] ) /* Min X */
589
newShape->param.poly.xmin = coords[i];
590
if( newShape->param.poly.xmax < coords[i] ) /* Max X */
591
newShape->param.poly.xmax = coords[i];
592
i++;
593
if( newShape->param.poly.ymin > coords[i] ) /* Min Y */
594
newShape->param.poly.ymin = coords[i];
595
if( newShape->param.poly.ymax < coords[i] ) /* Max Y */
596
newShape->param.poly.ymax = coords[i];
597
i++;
598
}
599
break;
600
}
601
602
} /* End of while( *currLoc ) */
603
/*
604
if (coords)printf("%.8f %.8f %.8f %.8f %.8f\n",
605
coords[0],coords[1],coords[2],coords[3],coords[4]);
606
*/
607
} /* End of if...else parse line */
608
} /* End of while( fgets(rgnFile) ) */
609
610
611
error:
612
613
if( *status )
614
fits_free_region( aRgn );
615
else
616
*Rgn = aRgn;
617
618
fclose( rgnFile );
619
free( currLine );
620
621
return( *status );
622
}
623
624
/*---------------------------------------------------------------------------*/
625
int fftrgn( double X,
626
double Y,
627
SAORegion *Rgn )
628
/* Test if the given point is within the region described by Rgn. X and */
629
/* Y are in pixel coordinates. */
630
/*---------------------------------------------------------------------------*/
631
{
632
double x, y, dx, dy, xprime, yprime, r;
633
RgnShape *Shapes;
634
int i;
635
int result = 0;
636
637
Shapes = Rgn->Shapes;
638
639
/* if an excluded region is given first, then implicitly */
640
/* assume a previous shape that includes the entire image. */
641
if (!Shapes->sign)
642
result = 1;
643
644
for( i=0; i<Rgn->nShapes; i++, Shapes++ ) {
645
646
/* only need to test if */
647
/* the point is not already included and this is an include region, */
648
/* or the point is included and this is an excluded region */
649
650
if ( (!result && Shapes->sign) || (result && !Shapes->sign) ) {
651
652
result = 1;
653
654
switch( Shapes->shape ) {
655
656
case box_rgn:
657
/* Shift origin to center of region */
658
xprime = X - Shapes->param.gen.p[0];
659
yprime = Y - Shapes->param.gen.p[1];
660
661
/* Rotate point to region's orientation */
662
x = xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
663
y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
664
665
dx = 0.5 * Shapes->param.gen.p[2];
666
dy = 0.5 * Shapes->param.gen.p[3];
667
if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
668
result = 0;
669
break;
670
671
case rectangle_rgn:
672
/* Shift origin to center of region */
673
xprime = X - Shapes->param.gen.p[5];
674
yprime = Y - Shapes->param.gen.p[6];
675
676
/* Rotate point to region's orientation */
677
x = xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
678
y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
679
680
dx = Shapes->param.gen.a;
681
dy = Shapes->param.gen.b;
682
if( (x < -dx) || (x > dx) || (y < -dy) || (y > dy) )
683
result = 0;
684
break;
685
686
case diamond_rgn:
687
/* Shift origin to center of region */
688
xprime = X - Shapes->param.gen.p[0];
689
yprime = Y - Shapes->param.gen.p[1];
690
691
/* Rotate point to region's orientation */
692
x = xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
693
y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
694
695
dx = 0.5 * Shapes->param.gen.p[2];
696
dy = 0.5 * Shapes->param.gen.p[3];
697
r = fabs(x/dx) + fabs(y/dy);
698
if( r > 1 )
699
result = 0;
700
break;
701
702
case circle_rgn:
703
/* Shift origin to center of region */
704
x = X - Shapes->param.gen.p[0];
705
y = Y - Shapes->param.gen.p[1];
706
707
r = x*x + y*y;
708
if ( r > Shapes->param.gen.a )
709
result = 0;
710
break;
711
712
case annulus_rgn:
713
/* Shift origin to center of region */
714
x = X - Shapes->param.gen.p[0];
715
y = Y - Shapes->param.gen.p[1];
716
717
r = x*x + y*y;
718
if ( r < Shapes->param.gen.a || r > Shapes->param.gen.b )
719
result = 0;
720
break;
721
722
case sector_rgn:
723
/* Shift origin to center of region */
724
x = X - Shapes->param.gen.p[0];
725
y = Y - Shapes->param.gen.p[1];
726
727
if( x || y ) {
728
r = atan2( y, x ) * 180.0 / myPI;
729
if( Shapes->param.gen.p[2] <= Shapes->param.gen.p[3] ) {
730
if( r < Shapes->param.gen.p[2] || r > Shapes->param.gen.p[3] )
731
result = 0;
732
} else {
733
if( r < Shapes->param.gen.p[2] && r > Shapes->param.gen.p[3] )
734
result = 0;
735
}
736
}
737
break;
738
739
case ellipse_rgn:
740
/* Shift origin to center of region */
741
xprime = X - Shapes->param.gen.p[0];
742
yprime = Y - Shapes->param.gen.p[1];
743
744
/* Rotate point to region's orientation */
745
x = xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
746
y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
747
748
x /= Shapes->param.gen.p[2];
749
y /= Shapes->param.gen.p[3];
750
r = x*x + y*y;
751
if( r>1.0 )
752
result = 0;
753
break;
754
755
case elliptannulus_rgn:
756
/* Shift origin to center of region */
757
xprime = X - Shapes->param.gen.p[0];
758
yprime = Y - Shapes->param.gen.p[1];
759
760
/* Rotate point to outer ellipse's orientation */
761
x = xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
762
y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
763
764
x /= Shapes->param.gen.p[4];
765
y /= Shapes->param.gen.p[5];
766
r = x*x + y*y;
767
if( r>1.0 )
768
result = 0;
769
else {
770
/* Repeat test for inner ellipse */
771
x = xprime * Shapes->param.gen.b + yprime * Shapes->param.gen.a;
772
y = -xprime * Shapes->param.gen.a + yprime * Shapes->param.gen.b;
773
774
x /= Shapes->param.gen.p[2];
775
y /= Shapes->param.gen.p[3];
776
r = x*x + y*y;
777
if( r<1.0 )
778
result = 0;
779
}
780
break;
781
782
case line_rgn:
783
/* Shift origin to first point of line */
784
xprime = X - Shapes->param.gen.p[0];
785
yprime = Y - Shapes->param.gen.p[1];
786
787
/* Rotate point to line's orientation */
788
x = xprime * Shapes->param.gen.cosT + yprime * Shapes->param.gen.sinT;
789
y = -xprime * Shapes->param.gen.sinT + yprime * Shapes->param.gen.cosT;
790
791
if( (y < -0.5) || (y >= 0.5) || (x < -0.5)
792
|| (x >= Shapes->param.gen.a) )
793
result = 0;
794
break;
795
796
case point_rgn:
797
/* Shift origin to center of region */
798
x = X - Shapes->param.gen.p[0];
799
y = Y - Shapes->param.gen.p[1];
800
801
if ( (x<-0.5) || (x>=0.5) || (y<-0.5) || (y>=0.5) )
802
result = 0;
803
break;
804
805
case poly_rgn:
806
if( X<Shapes->param.poly.xmin || X>Shapes->param.poly.xmax
807
|| Y<Shapes->param.poly.ymin || Y>Shapes->param.poly.ymax )
808
result = 0;
809
else
810
result = Pt_in_Poly( X, Y, Shapes->param.poly.nPts,
811
Shapes->param.poly.Pts );
812
break;
813
}
814
815
if( !Shapes->sign ) result = !result;
816
817
}
818
}
819
820
return( result );
821
}
822
823
/*---------------------------------------------------------------------------*/
824
void fffrgn( SAORegion *Rgn )
825
/* Free up memory allocated to hold the region data. */
826
/*---------------------------------------------------------------------------*/
827
{
828
int i;
829
830
for( i=0; i<Rgn->nShapes; i++ )
831
if( Rgn->Shapes[i].shape == poly_rgn )
832
free( Rgn->Shapes[i].param.poly.Pts );
833
if( Rgn->Shapes )
834
free( Rgn->Shapes );
835
free( Rgn );
836
}
837
838
/*---------------------------------------------------------------------------*/
839
static int Pt_in_Poly( double x,
840
double y,
841
int nPts,
842
double *Pts )
843
/* Internal routine for testing whether the coordinate x,y is within the */
844
/* polygon region traced out by the array Pts. */
845
/*---------------------------------------------------------------------------*/
846
{
847
int i, j, flag=0;
848
double prevX, prevY;
849
double nextX, nextY;
850
double dx, dy, Dy;
851
852
nextX = Pts[nPts-2];
853
nextY = Pts[nPts-1];
854
855
for( i=0; i<nPts; i+=2 ) {
856
prevX = nextX;
857
prevY = nextY;
858
859
nextX = Pts[i];
860
nextY = Pts[i+1];
861
862
if( (y>prevY && y>=nextY) || (y<prevY && y<=nextY)
863
|| (x>prevX && x>=nextX) )
864
continue;
865
866
/* Check to see if x,y lies right on the segment */
867
868
if( x>=prevX || x>nextX ) {
869
dy = y - prevY;
870
Dy = nextY - prevY;
871
872
if( fabs(Dy)<1e-10 ) {
873
if( fabs(dy)<1e-10 )
874
return( 1 );
875
else
876
continue;
877
}
878
879
dx = prevX + ( (nextX-prevX)/(Dy) ) * dy - x;
880
if( dx < -1e-10 )
881
continue;
882
if( dx < 1e-10 )
883
return( 1 );
884
}
885
886
/* There is an intersection! Make sure it isn't a V point. */
887
888
if( y != prevY ) {
889
flag = 1 - flag;
890
} else {
891
j = i+1; /* Point to Y component */
892
do {
893
if( j>1 )
894
j -= 2;
895
else
896
j = nPts-1;
897
} while( y == Pts[j] );
898
899
if( (nextY-y)*(y-Pts[j]) > 0 )
900
flag = 1-flag;
901
}
902
903
}
904
return( flag );
905
}
906
Loggerhead is a web-based interface for Breezy
Version: 2.0.1