diff options
Diffstat (limited to 'funtools/filter/imregions_c.tmpl')
| -rw-r--r-- | funtools/filter/imregions_c.tmpl | 2585 |
1 files changed, 2585 insertions, 0 deletions
diff --git a/funtools/filter/imregions_c.tmpl b/funtools/filter/imregions_c.tmpl new file mode 100644 index 0000000..639fefc --- /dev/null +++ b/funtools/filter/imregions_c.tmpl @@ -0,0 +1,2585 @@ +/* + NB: MAKE SURE YOU EDIT THE TEMPLATE FILE!!!! +*/ + +#ifndef FILTER_PTYPE +#include <regions.h> +#endif + +#ifndef UNUSED +#ifdef __GNUC__ +# define UNUSED(x) UNUSED_ ## x __attribute__((__unused__)) +#else +# define UNUSED(x) UNUSED_ ## x +#endif +#endif + +/* panda and pie incorrectly used astronomical angles. fixed 4/2004 */ +#define USE_ASTRO_ANGLE 0 + +/* we might want to avoid (x86) extended precision problems */ +#define USE_FPU_DOUBLE _FPUDBL_ +#if USE_FPU_DOUBLE +#include <fpu_control.h> +static fpu_control_t _cw; +#define FPU_DOUBLE {fpu_control_t _cw2; _FPU_GETCW(_cw); _cw2 = _cw & ~_FPU_EXTENDED; _cw2 |= _FPU_DOUBLE; _FPU_SETCW(_cw2);} +#define FPU_RESTORE {_FPU_SETCW(_cw);} +#else +#define FPU_DOUBLE +#define FPU_RESTORE +#endif + +/* add to FilterOpen to cause this module to be loaded for dynamic linking */ +static int imregno=0; +void initimregions(void) +{ + imregno++; + return; +} + +static void markx(GFilt g, int UNUSED(sno), int flag, int type, int x, int y) +{ + /* don't mark exclude regions */ + if( type == TOK_EREG ) + return; + /* for include, we might extend the x limits */ + if( flag ){ + if( x <= g->x0s[y] ) + g->x0s[y] = max(x,g->x0); + if( x >= g->x1s[y] ) + g->x1s[y] = min(x,g->x1); + } + /* for exclude, we have to look at the full line */ + else{ + g->x0s[y] = g->x0; + g->x1s[y] = g->x1; + } +} + +static void marky(GFilt g, int sno, int flag, int type) +{ + int i; + /* don't mark exclude regions */ + if( type == TOK_EREG ) + return; + /* mark include shape */ + if( flag ){ + for(i=max(g->y0,g->shapes[sno].ystart); + i<=min(g->y1,g->shapes[sno].ystop); + i++) + g->ybuf[i] = 1; + } + /* mark exclude shape */ + else{ + for(i=g->y0; i<=g->shapes[sno].ystart-1; i++) + g->ybuf[i] = 1; + for(i=g->shapes[sno].ystop+1; i<=g->y1; i++) + g->ybuf[i] = 1; + } +} + +static int +imagemaskcmp(const void *s1, const void *s2) +{ + FilterMask f1 = (FilterMask)s1; + FilterMask f2 = (FilterMask)s2; + + if( f1->y < f2->y ){ + return -1; + } + else if( f1->y > f2->y ){ + return 1; + } + else{ + if( f1->xstart < f2->xstart ){ + return -1; + } + else{ + return 1; + } + } +} + +/* ***************************** shape support ***************************** */ + +static void quadeq(double a, double b, double c, + double *x1, double *x2, int *nr, int *nc) +{ + double dis, q; + if( feq(a,0.0) ){ + *nc = 0; + if( feq(b,0.0) ){ + *nr = 0; *x1 = 0.0; + } + else{ + *nr = 1; *x1 = -c / b; + } + *x2 = *x1; + } + else{ + dis = b*b - 4.0 * a * c; + if( dis > 0.0 ){ + *nr = 2; *nc = 0; + dis = sqrt(dis); + if( b < 0.0 ) dis = -dis; + q = -0.5 * (b + dis); + *x1 = q/a; *x2 = c/q; + if(*x1 > *x2){ + q = *x1; *x1 = *x2; *x2 = q; + } + } + else if( feq(dis,0.0) ){ + *nr = 1; *nc = 0; *x1 = - 0.5 * b / a; *x2 = *x1; + } + else{ + *nr = 0; *nc = 2; *x1 = - 0.5 * b / a; *x2 = 0.5 * sqrt(-dis) / a; + } + } +} + +static void rgs_mark(GFilt g, Scan *scanlist, + int sno, int flag, int type, int xval, int yval ) +{ + Scan scanmark, mark; + /* since yval is used as an index, make sure its within limits */ + if(yval < g->y0) yval = g->y0; + if(yval > g->y1) yval = g->y1; + mark = (Scan)calloc(1, sizeof(ScanRec)); + mark->x = xval; + /* sanity check */ + if( !scanlist ) return; + /* starts are installed at back of list for given x */ + if( !scanlist[yval] || ((scanlist[yval])->x > xval) ){ + mark->next = scanlist[yval]; + scanlist[yval] = mark; + } else { + scanmark = scanlist[yval]; + while( scanmark->next && ((scanmark->next)->x < xval) ) + scanmark = scanmark->next; + mark->next = scanmark->next; + scanmark->next = mark; + } + markx(g, sno, flag, type, xval, yval); +} + +static void rgs_segment(GFilt g, Scan *scanlist, int width, int height, + int sno, int flag, int type, + double x1, double y1, double x2, double y2) +{ + int ystart, ystop, yval, xval; + double invslope, xoffset; + + ystart = PIXINCL(y1); + if( ystart < 1 ) ystart = 1; + /* note: PIXINCL(stop) is 1st pixel not counted */ + ystop = PIXINCL(y2) - 1; + if( ystop > height ) ystop = height; + /* ignore segment if there is no positive slope in integer coords */ + if( (ystart > ystop) || (ystop < 1) ) + return; + /* use inverse slope (run/rise) to get x given y with a multiply */ + invslope = (x1 - x2) / (y1 - y2); + xoffset = x1 + ((ystart - y1) * invslope); + for(yval=ystart; yval<=ystop; yval++){ + xval = PIXINCL(xoffset); + /* clip line to edges of image area (actually bend line) */ + if(xval < 1) xval = 1; + if(xval > width) xval = width + 1; + rgs_mark(g, scanlist, sno, flag, type, xval, yval); + xoffset = xoffset + invslope; + } +} + +static void _polygoni(GFilt g, int qt, int UNUSED(rno), int sno, int flag, + int type, double *vx, double *vy, int count) +{ + int i, j; + double xlo, xhi; + double ylo, yhi; + +#ifdef ALREADY_DONE + /* divide by block factor to get real endpoints */ + for(i=0; i<count; i++){ + vx[i] = (vx[i] - g->xmin)/g->block + 1.0; + vy[i] = (vy[i] - g->ymin)/g->block + 1.0; + } +#endif + /* find the limits */ + xlo = vx[0]; + xhi = xlo; + ylo = vy[0]; + yhi = ylo; + for(i=0; i<count; i++){ + if(vx[i] > xhi) xhi = vx[i]; + if(vx[i] < xlo) xlo = vx[i]; + if(vy[i] > yhi) yhi = vy[i]; + if(vy[i] < ylo) ylo = vy[i]; + } + if( qt && (sno > 1) ){ + g->shapes[sno].ystart = g->shapes[sno-1].ystart; + g->shapes[sno].ystop = g->shapes[sno-1].ystop; + } + else{ + g->shapes[sno].ystart = max(g->y0,PIXINCL(ylo)); + g->shapes[sno].ystop = min(g->y1,PIXINCL(yhi) - 1); + } + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + marky(g, sno, flag, type); + + /* mark all horizontal segment crossings */ + /* start with segment between last and first point */ + j = count-1; + for(i=0; i<count; i++){ + /* make segments always run from lower y to higher y */ + if(vy[i] > vy[j]){ + rgs_segment(g, g->shapes[sno].scanlist, g->x1, g->y1, + sno, flag, type, vx[j], vy[j], vx[i], vy[i]); + } + else{ + rgs_segment(g, g->shapes[sno].scanlist, g->x1, g->y1, + sno, flag, type, vx[i], vy[i], vx[j], vy[j]); + } + j = i; + } +} + +static int corner_vertex(int index, int width, int height, + double *x, double *y) +{ + switch (index) { + case 1: + *x = 0.0; + *y = height + 1; + break; + case 2: + *x = 0.0; + *y = 0.0; + break; + case 3: + *x = width + 1; + *y = 0.0; + break; + case 4: + *x = width + 1; + *y = height + 1; + default: + break; + } + index = index + 1; + if(index > 4) index = 1; + return(index); +} + +static int pie_intercept(int width, int height, double xcen, double ycen, + double angle, double *xcept, double *ycept) +{ + double angl, slope; /* l: angle and slope of ray */ + angl = angle; + /* put angles in normal range */ + while (angl < 0.0) + angl = angl + 360.0; + while (angl >= 360.0) + angl = angl - 360.0; + /* check for a horizontal angle */ +#if USE_ASTRO_ANGLE + if(fabs(angl - 90.0) < SMALL_NUMBER) { +#else + if(fabs(angl - 180.0) < SMALL_NUMBER) { +#endif + *xcept = 0.0; + *ycept = ycen; + return(2); + } +#if USE_ASTRO_ANGLE + if(fabs(angl - 270.0) < SMALL_NUMBER) { +#else + if(fabs(angl - 0.0) < SMALL_NUMBER) { +#endif + *xcept = width + 1; + *ycept = ycen; + return(4); + } +#if USE_ASTRO_ANGLE + /* convert to a Cartesian angle */ + angl = angl + 90.0; +#endif + if(angl >= 360.0) + angl = angl - 360.0; + if(angl < 180.0) { + *ycept = height + 1; + /* rule out vertical line */ + if(fabs(angl - 90.0) < SMALL_NUMBER) { + *xcept = xcen; + return(1); + } + } else { + *ycept = 0.0; + /* rule out vertical line */ + if(fabs(angl - 270.0) < SMALL_NUMBER) { + *xcept = xcen; + return(3); + } + } + /* convert to radians */ + angl = (angl / 180.0) * M_PI; + /* calculate slope */ + slope = tan(angl); + /* calculate intercept with designated y edge */ + *xcept = xcen + ((*ycept - ycen) / slope); + if(*xcept < 0) { + *ycept = (ycen - (xcen * slope)); + *xcept = 0.0; + return(2); + } else if(*xcept > (width + 1)) { + *ycept = (ycen + ((width + 1 - xcen) * slope)); + *xcept = width + 1; + return(4); + } else { + if(*ycept < height) + return(3); + else + return(1); + } +} + +void _impiei(GFilt g, int qt, int rno, int sno, int flag, int type, + double UNUSED(x), double UNUSED(y), + double xcen, double ycen, double angle1, double angle2) +{ + int width, height; /* l: image mask width and height */ + double sweep; /* l: sweep between cut angles */ + double vx[7], vy[7]; /* l: arrays of vertices for polygon */ + int count; /* l: number of polygon vertices */ + int intrcpt1, intrcpt2; /* l: side intercepted by each cut */ + double x2, y2; /* l: coordinates of second intercept */ + + /* NB: do not use x and y variables, they have bogus values */ + /* divide by block factor to get "real" parameters */ + xcen = (xcen - g->xmin)/g->block + 1.0; + ycen = (ycen - g->ymin)/g->block + 1.0; + /* temps */ + width = g->x1; + height = g->y1; + /* start listing vertices of polygon */ + vx[0] = xcen; + vy[0] = ycen; + sweep = angle2 - angle1; + /* if sweep is too small to be noticed, don't bother */ + if(fabs(sweep) < SMALL_NUMBER) + return; + if (sweep < 0.0) sweep = sweep + 360.0; + intrcpt1 = pie_intercept(width, height, xcen, ycen, angle1, + &(vx[1]), &(vy[1])); + intrcpt2 = pie_intercept(width, height, xcen, ycen, angle2, + &x2, &y2); + count = 2; + /* if angles intercept same side and slice is between them, no corners */ + /* else, mark corners until reaching side with second angle intercept */ + if((intrcpt1 != intrcpt2) || (sweep > 180.0)){ + do{ + intrcpt1 = corner_vertex(intrcpt1, width, height, + &(vx[count]), &(vy[count])); + count = count + 1; + }while(intrcpt1 != intrcpt2); + } + vx[count] = x2; + vy[count] = y2; + count++; + _polygoni(g, qt, rno, sno, flag, type, vx, vy, count); +} + +/* ***************************** shapes ********************************** */ + +void imannulusi(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, double iradius, double oradius) +{ + int yy; + double dval; + double xoff, yoff; + Scan *scanlist; + + /* NB: do not use x and y variables, they have bogus values */ + + if( iradius == 0 ){ + imcirclei(g, rno, sno, flag, type, x, y, xcen, ycen, oradius); + return; + } + /* divide by block factor to get "real" parameters */ + xcen = (xcen - g->xmin)/g->block + 1.0; + ycen = (ycen - g->ymin)/g->block + 1.0; + iradius /= (double)g->block; + oradius /= (double)g->block; + /* set y limits */ + if( PIXSTART(ycen - oradius) < PIXSTOP(ycen + oradius) ){ + g->shapes[sno].ystart = max(g->y0,PIXSTART(ycen - oradius)); + g->shapes[sno].ystop = min(g->y1,PIXSTOP(ycen + oradius)); + } + else{ + g->shapes[sno].ystart = min(g->y1,PIXSTOP(ycen + oradius)); + g->shapes[sno].ystop = max(g->y0,PIXSTART(ycen - oradius)); + } + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + scanlist = g->shapes[sno].scanlist; + marky(g, sno, flag, type); + /* calculate start/stop values for each y line */ + for(yy=g->shapes[sno].ystart; yy<=g->shapes[sno].ystop; yy++){ + yoff = PIXCEN(yy) - ycen; + if( (dval=(oradius * oradius) - (yoff * yoff)) > 0.0 ){ + xoff = sqrt(dval); + if( PIXSTART(xcen - xoff) <= PIXSTOP(xcen + xoff) ){ + rgs_mark(g, scanlist, sno, flag, type, PIXSTART(xcen - xoff), yy); + rgs_mark(g, scanlist, sno, flag, type, PIXSTOP(xcen + xoff), yy); + } + if( (dval=(iradius * iradius) - (yoff * yoff)) > 0.0 ){ + xoff = sqrt(dval); + if( PIXSTART(xcen - xoff) <= PIXSTOP(xcen + xoff) ){ + rgs_mark(g, scanlist, sno, flag, type, PIXSTART(xcen - xoff), yy); + rgs_mark(g, scanlist, sno, flag, type, PIXSTOP(xcen + xoff), yy); + } + } + } + } +} + +int imannulus(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, double iradius, double oradius) +{ + Scan scan; + + if( iradius == 0 ){ + return(imcircle(g, rno, sno, flag, type, x, y, xcen, ycen, oradius)); + } + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + if( flag ){ + if(y < g->shapes[sno].ystart ) return 0; + if(y > g->shapes[sno].ystop ) return 0; + } + scan = g->shapes[sno].scanlist[(int)y]; + if( (scan && + ((y >= g->shapes[sno].ystart) && (y <= g->shapes[sno].ystop)) && + (scan->next->next ? + (((x >= scan->x) && (x <= scan->next->next->next->x)) && + !((x >= scan->next->x) && (x <= scan->next->next->x))) : + ((x >= scan->x) && (x <= scan->next->x)))) + ==flag ){ + if( rno && flag ) g->rid = rno; + return 1; + } + else + return 0; +} + +void imboxi(GFilt g, int rno, int sno, int flag, int type, + double UNUSED(x), double UNUSED(y), + double xcen, double ycen, double xwidth, double yheight, + double angle) +{ + double angl; /* l: Cartesian angle in radians */ + double half_width, half_height;/* l: radii (1/2 width and height) */ + double cosangl, sinangl; /* l: sine, cosine of the Cartesian angle */ + double hw_cos, hw_sin; /* l: products of half_width with sin, cos */ + double hh_cos, hh_sin; /* l: products of half_height with sin, cos */ + double cornerx[4], cornery[4]; /* l: arrays of x and y coords of 4 corners */ + + /* NB: do not use x and y variables, they have bogus values */ + if( (xwidth == 0) && (yheight==0) ){ + return; + } + /* divide by block factor to get "real" parameters */ + xcen = (xcen - g->xmin)/g->block + 1.0; + ycen = (ycen - g->ymin)/g->block + 1.0; + xwidth /= (double)g->block; + yheight /= (double)g->block; +#if USE_ASTRO_ANGLE + /* convert to a Cartesian angle; save angle for use in multi or slices */ + angl = angle + 90.0; +#else + angl = angle; +#endif + while (angl >= 360.0) angl = angl - 360.0; + /* convert to radians */ + angl = (angl / 180.0) * M_PI; + sinangl = sin (angl); + cosangl = cos (angl); +#if USE_ASTRO_ANGLE + /* since we rotate by 90.0 degrees to get from astro angle to cartesian, */ + /* we also need to switch the width and height. we do this secretly so */ + /* that the display will turn out right, by doing it in the half terms */ + half_width = yheight / 2.0; + half_height = xwidth / 2.0; +#else + half_width = xwidth / 2.0; + half_height = yheight / 2.0; +#endif + hw_cos = half_width * cosangl; + hw_sin = half_width * sinangl; + hh_cos = half_height * cosangl; + hh_sin = half_height * sinangl; +#if USE_ASTRO_ANGLE + cornerx[0] = xcen - hw_cos - hh_sin; + cornery[0] = ycen - hw_sin + hh_cos; + cornerx[1] = xcen + hw_cos - hh_sin; + cornery[1] = ycen + hw_sin + hh_cos; + cornerx[2] = xcen + hw_cos + hh_sin; + cornery[2] = ycen + hw_sin - hh_cos; + cornerx[3] = xcen - hw_cos + hh_sin; + cornery[3] = ycen - hw_sin - hh_cos; +#else + cornerx[0] = xcen - hw_cos + hh_sin; + cornery[0] = ycen - hh_cos - hw_sin; + cornerx[1] = xcen - hw_cos - hh_sin; + cornery[1] = ycen + hh_cos - hw_sin; + cornerx[2] = xcen + hw_cos - hh_sin; + cornery[2] = ycen + hh_cos + hw_sin; + cornerx[3] = xcen + hw_cos + hh_sin; + cornery[3] = ycen - hh_cos + hw_sin; +#endif + _polygoni(g, 0, rno, sno, flag, type, cornerx, cornery, 4); +} + +int imbox(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double UNUSED(xcen), double UNUSED(ycen), + double xwidth, double yheight, + double UNUSED(angle)) +{ + if( (xwidth == 0) && (yheight==0) ){ + return(!flag); + } + return impolygon(g, rno, sno, flag, type, x, y); +} + +void imcirclei(GFilt g, int UNUSED(rno), int sno, int flag, int type, + double UNUSED(x), double UNUSED(y), + double xcen, double ycen, double radius) +{ + int yy; + double dval; + double xoff, yoff; + Scan *scanlist; + + /* NB: do not use x and y variables, they have bogus values */ + if( radius == 0 ){ + return; + } + /* divide by block factor to get "real" parameters */ + xcen = (xcen - g->xmin)/g->block + 1.0; + ycen = (ycen - g->ymin)/g->block + 1.0; + radius /= (double)g->block; + /* set y limits */ + if( PIXSTART(ycen - radius) < PIXSTOP(ycen + radius) ){ + g->shapes[sno].ystart = max(g->y0,PIXSTART(ycen - radius)); + g->shapes[sno].ystop = min(g->y1,PIXSTOP(ycen + radius)); + } + else{ + g->shapes[sno].ystart = min(g->y1,PIXSTOP(ycen + radius)); + g->shapes[sno].ystop = max(g->y0,PIXSTART(ycen - radius)); + } + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + scanlist = g->shapes[sno].scanlist; + marky(g, sno, flag, type); + /* calculate start/stop values for each y line */ + for(yy=g->shapes[sno].ystart; yy<=g->shapes[sno].ystop; yy++){ + yoff = PIXCEN(yy) - ycen; + if( (dval=(radius * radius) - (yoff * yoff))>=0.0 ){ + xoff = sqrt(dval); + if( PIXSTART(xcen - xoff) <= PIXSTOP(xcen + xoff) ){ + rgs_mark(g, scanlist, sno, flag, type, PIXSTART(xcen - xoff), yy); + rgs_mark(g, scanlist, sno, flag, type, PIXSTOP(xcen + xoff), yy); + } + } + } +} + +int imcircle(GFilt g, int rno, int sno, int flag, int UNUSED(type), + double x, double y, + double UNUSED(xcen), double UNUSED(ycen), double radius) +{ + Scan scan; + + if( radius == 0 ){ + return(!flag); + } + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + if( flag ){ + if(y < g->shapes[sno].ystart ) return 0; + if(y > g->shapes[sno].ystop ) return 0; + } + scan = g->shapes[sno].scanlist[(int)y]; + if( (scan && + ((y>=g->shapes[sno].ystart) && (y<=g->shapes[sno].ystop)) && + ((x >= scan->x) && (x <= (scan->next)->x))) == flag ){ + if( rno && flag ) g->rid = rno; + return 1; + } + else{ + return 0; + } +} + +void imellipsei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, double xrad, double yrad, + double angle) +{ + int yy; + int nr, nc; + double yhi; + double yoff; + double xboff, xfoff; + double angl; + double sinangl, cosangl; + double cossq, sinsq, xradsq, yradsq; + double a, b_partial, c_partial; + double b, c; + Scan *scanlist; + + /* NB: do not use x and y variables, they have bogus values */ + + if( xrad == yrad ){ + imcirclei(g, rno, sno, flag, type, x, y, xcen, ycen, xrad); + return; + } + /* divide by block factor to get "real" parameters */ + xcen = (xcen - g->xmin)/g->block + 1.0; + ycen = (ycen - g->ymin)/g->block + 1.0; + xrad /= (double)g->block; + yrad /= (double)g->block; + /* set worst case limits (xrad axis parallel to vertical axis) */ +#if USE_ASTRO_ANGLE + /* convert to a Cartesian angle; save "angle" for use by other routines */ + angl = angle + 90.0; +#else + angl = angle; +#endif + while( angl >= 360.0 ) + angl = angl - 360.0; + /* convert to radians */ + angl = (angl / 180.0) * M_PI; + sinangl = sin(angl); + cosangl = cos(angl); + /* calculate approximate y limits */ + /* choose lesser of containing rotbox and circle */ + FPU_DOUBLE +#if USE_ASTRO_ANGLE + yhi = fabs(sinangl * yrad) + fabs(cosangl * xrad); +#else + yhi = fabs(sinangl * xrad) + fabs(cosangl * yrad); +#endif + yhi = min(yhi, max(yrad, xrad)); + /* set y limits */ + if( PIXSTART(ycen - yhi) < PIXSTOP(ycen + yhi) ){ + g->shapes[sno].ystart = max(g->y0,PIXSTART(ycen - yhi)); + g->shapes[sno].ystop = min(g->y1,PIXSTOP(ycen + yhi)); + } + else{ + g->shapes[sno].ystart = min(g->y1,PIXSTOP(ycen + yhi)); + g->shapes[sno].ystop = max(g->y0,PIXSTART(ycen - yhi)); + } + FPU_RESTORE + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + scanlist = g->shapes[sno].scanlist; + marky(g, sno, flag, type); + /* prepare partials for quadratic equation solutions to coordinates */ + cossq = cosangl * cosangl; + sinsq = sinangl * sinangl; +#if USE_ASTRO_ANGLE + /* because we rotate by 90.0 degrees to get from astro angle to */ + /* cartesian, we also need to switch the x and y axes. we do this */ + /* secretly so that the display will turn out right, by doing it in */ + /* the sq terms */ + xradsq = yrad * yrad; + yradsq = xrad * xrad; +#else + xradsq = xrad * xrad; + yradsq = yrad * yrad; +#endif + /* fill in as much of a,b,c as we can */ + a = (cossq / xradsq) + (sinsq / yradsq); + b_partial = (2.0 * sinangl) * ((cosangl / xradsq) - (cosangl / yradsq)); + c_partial = (sinsq / xradsq) + (cossq / yradsq); + /* calculate start/stop values for each y line */ + for(yy=g->shapes[sno].ystart; yy<=g->shapes[sno].ystop; yy++){ + yoff = yy - ycen; + b = b_partial * yoff; + c = (c_partial * yoff * yoff) - 1.0; + /* solve quadratic */ + quadeq (a, b, c, &xboff, &xfoff, &nr, &nc); + /* if real roots */ + if( nr != 0 ) { + /* translate x coordinates */ + rgs_mark(g, scanlist, sno, flag, type, PIXSTART(xcen + xboff), yy); + rgs_mark(g, scanlist, sno, flag, type, PIXSTOP(xcen + xfoff), yy); + } + } +} + +int imellipse(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, double xrad, double yrad, + double UNUSED(angle)) +{ + Scan scan; + + if( xrad == yrad ){ + return(imcircle(g, rno, sno, flag, type, x, y, xcen, ycen, xrad)); + } + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + if( flag ){ + if(y < g->shapes[sno].ystart ) return 0; + if(y > g->shapes[sno].ystop ) return 0; + } + scan = g->shapes[sno].scanlist[(int)y]; + FPU_DOUBLE + if( (scan && + ((y>=g->shapes[sno].ystart) && (y<=g->shapes[sno].ystop)) && + ((x >= scan->x) && (x <= (scan->next)->x))) == flag ){ + if( rno && flag ) g->rid = rno; + FPU_RESTORE + return 1; + } + else{ + FPU_RESTORE + return 0; + } +} + +void imfieldi(GFilt g, int UNUSED(rno), int sno, int flag, int type, + double UNUSED(x), double UNUSED(y)) +{ + int yy; + Scan *scanlist; + + /* NB: do not use x and y variables, they have bogus values */ + g->shapes[sno].ystart = g->y0; + g->shapes[sno].ystop = g->y1; + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + scanlist = g->shapes[sno].scanlist; + marky(g, sno, flag, type); + /* calculate start/stop values for each y line */ + for(yy=g->shapes[sno].ystart; yy<=g->shapes[sno].ystop; yy++){ + rgs_mark(g, scanlist, sno, flag, type, g->x0, yy); + rgs_mark(g, scanlist, sno, flag, type, g->x1, yy); + } +} + +int imfield(GFilt g, int rno, int UNUSED(sno), int flag, int UNUSED(type), + double UNUSED(x), double UNUSED(y)) +{ + if( flag ){ + if( rno && flag ) g->rid = rno; + return 1; + } + else{ + return 0; + } +} + +void imlinei(GFilt g, int UNUSED(rno), int sno, int flag, int type, + double UNUSED(x), double UNUSED(y), + double x1, double y1, double x2, double y2) +{ + double vx[2]; + double vy[2]; + int xval, yval; + double invslope, xoffset; + + /* NB: do not use x and y variables, they have bogus values */ + /* divide by block factor to get "real" parameters */ + x1 = (x1 - g->xmin)/g->block + 1.0; + y1 = (y1 - g->ymin)/g->block + 1.0; + x2 = (x2 - g->xmin)/g->block + 1.0; + y2 = (y2 - g->ymin)/g->block + 1.0; + /* order by increasing y */ + if( y1 < y2 ){ + vx[0] = x1; + vy[0] = y1; + vx[1] = x2; + vy[1] = y2; + } + else{ + vx[0] = x2; + vy[0] = y2; + vx[1] = x1; + vy[1] = y1; + } + /* set y limits */ + g->shapes[sno].ystart = PIXNUM(vy[0]); + g->shapes[sno].ystop = PIXNUM(vy[1]); + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + marky(g, sno, flag, type); + if( feq(vy[0],vy[1]) ){ + rgs_mark(g, g->shapes[sno].scanlist, sno, flag, type, + (int)min(vx[0],vx[1]), (int)vy[0]); + rgs_mark(g, g->shapes[sno].scanlist, sno, flag, type, + (int)max(vx[0],vx[1]), (int)vy[0]); + } + else{ + /* mark all horizontal segment crossings */ + invslope = (vx[0] - vx[1]) / (vy[0] - vy[1]); + xoffset = vx[0]; + for(yval=vy[0]; yval<=vy[1]; yval++){ + xval = xoffset; + rgs_mark(g, g->shapes[sno].scanlist, sno, flag, type, xval, yval); + xoffset = xoffset + invslope; + } + } +} + +int imline(GFilt g, int rno, int sno, int flag, int UNUSED(type), + double x, double y, + double UNUSED(x1), double UNUSED(y1), + double UNUSED(x2), double UNUSED(y2)) +{ + Scan scan; + + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + if( flag ){ + if(y < g->shapes[sno].ystart ) return 0; + if(y > g->shapes[sno].ystop ) return 0; + } + scan = g->shapes[sno].scanlist[(int)y]; + if( (scan && + ((x==(int)scan->x) || + (scan->next&&((x>=(int)scan->x)&&(x<=(int)scan->next->x)))))==flag ) { + if( rno && flag ) g->rid = rno; + return 1; + } + else + return 0; +} + +void impiei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, double angle1, double angle2) +{ + _impiei(g, 0, rno, sno, flag, type, x, y, xcen, ycen, angle1, angle2); +} + +int impie(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double UNUSED(xcen), double UNUSED(ycen), + double UNUSED(angle1), double UNUSED(angle2)) +{ + return impolygon(g, rno, sno, flag, type, x, y); +} + +void imqtpiei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, double angle1, double angle2) +{ + _impiei(g, 1, rno, sno, flag, type, x, y, xcen, ycen, angle1, angle2); +} + +int imqtpie(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double UNUSED(xcen), double UNUSED(ycen), + double UNUSED(angle1), double UNUSED(angle2)) +{ + return impolygon(g, rno, sno, flag, type, x, y); +} + +void impointi(GFilt g, int UNUSED(rno), int sno, int flag, int type, + double UNUSED(x), double UNUSED(y), + double xcen, double ycen) +{ + /* NB: do not use x and y variables, they have bogus values */ + /* divide by block factor to get "real" parameters */ + xcen = (xcen - g->xmin)/g->block + 1.0; + ycen = (ycen - g->ymin)/g->block + 1.0; + /* set y limits */ + g->shapes[sno].ystart = PIXNUM(ycen); + g->shapes[sno].ystop = PIXNUM(ycen); + g->shapes[sno].scanlist = (Scan *)calloc(g->y1+1, sizeof(Scan)); + marky(g, sno, flag, type); + rgs_mark(g, g->shapes[sno].scanlist, sno, flag, type, + PIXNUM(xcen), PIXNUM(ycen)); +} + +int impoint(GFilt g, int rno, int sno, int flag, int UNUSED(type), + double x, double y, + double UNUSED(xcen), double UNUSED(ycen)) +{ + Scan scan; + + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + if( flag ){ + if(y < g->shapes[sno].ystart ) return 0; + if(y > g->shapes[sno].ystop ) return 0; + } + scan = g->shapes[sno].scanlist[(int)y]; + if( (scan && + (y == (int)g->shapes[sno].ystart) && + (x == (int)scan->x))==flag ) { + if( rno && flag ) g->rid = rno; + return 1; + } + else + return 0; +} + +#ifdef __STDC__ +void +impolygoni(GFilt g, int rno, int sno, int flag, int type, + double UNUSED(x), double y, ...) +{ + double *vx=NULL, *vy=NULL; + int count, maxcount; + va_list args; + va_start(args, y); +#else +void impolygoni(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double *vx=NULL, *vy=NULL; + int count, maxcount; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); +#endif + /* NB: do not use x and y variables, they have bogus values */ + /* allocate space for x,y arguments */ + maxcount = MASKINC; + vx = (double *)calloc(maxcount, sizeof(double)); + vy = (double *)calloc(maxcount, sizeof(double)); + /* look for x,y arguments */ + count = 0; + while( 1 ){ + if( (count + 1) >= maxcount ){ + maxcount += MASKINC; + vx = (double *)realloc(vx, maxcount*sizeof(double)); + vy = (double *)realloc(vy, maxcount*sizeof(double)); + } + vx[count] = va_arg(args, double); + vy[count] = va_arg(args, double); + if( feq(vx[count],PSTOP) && feq(vy[count],PSTOP) ) + break; + vx[count] = (vx[count] - g->xmin)/g->block + 1.0; + vy[count] = (vy[count] - g->ymin)/g->block + 1.0; + count ++; + } + va_end(args); + /* realloc to actual size */ + vx = (double *)realloc(vx, count*sizeof(double)); + vy = (double *)realloc(vy, count*sizeof(double)); + /* call the common routine */ + _polygoni(g, 0, rno, sno, flag, type, vx, vy, count); + if( vx ) free(vx); + if( vy ) free(vy); +} + +#ifdef __STDC__ +int +impolygon(GFilt g, int rno, int sno, int flag, int UNUSED(type), + double x, double y, ...) +{ + int crossings; + Scan scan; + va_list args; + va_start(args, y); +#else +int impolygon(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + int crossings; + Scan scan; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); +#endif + va_end(args); + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + if( flag ){ + if(y < g->shapes[sno].ystart ) return 0; + if(y > g->shapes[sno].ystop ) return 0; + } + /* no initialization of x for this row, just jump right in */ + if( (y>=g->shapes[sno].ystart)&&(y<=g->shapes[sno].ystop) ){ + crossings = 0; + for(scan=g->shapes[sno].scanlist[(int)y]; scan; scan=scan->next){ + if( x >= scan->x ) + crossings++; + else + break; + } + if( (crossings%2) == flag ){ + if( rno && flag ) g->rid = rno; + return 1; + } + else{ + return 0; + } + } + else{ + return !flag; + } +} + +void imnannulusi(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double lo, double hi, int n) +{ + int i; + int xsno; + double dinc; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + dinc = (hi - lo)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init all shapes */ + imannulusi(g, 0, xsno, flag, type, x, y, xcen, ycen, lo, hi); + for(i=0; i<n; i++){ + imannulusi(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lo+(i*dinc), lo+((i+1)*dinc)); + } +} + +void imnboxi(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double lox, double loy, double hix, double hiy, int n, + double ang) +{ + int i; + int xsno; + double dincx; + double dincy; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + dincx = (hix - lox)/n; + dincy = (hiy - loy)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init all shapes */ + imboxi(g, 0, xsno, flag, type, x, y, xcen, ycen, hix, hiy, ang); + imboxi(g, 0, xsno+1, flag, type, x, y, xcen, ycen, lox, loy, ang); + for(i=0; i<n; i++){ + imboxi(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lox+((i+1)*dincx), loy+((i+1)*dincy), ang); + } +} + +void imnellipsei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double lox, double loy, double hix, double hiy, int n, + double ang) +{ + int i; + int xsno; + double dincx; + double dincy; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + dincx = (hix - lox)/n; + dincy = (hiy - loy)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init all shapes */ + imellipsei(g, 0, xsno, flag, type, x, y, xcen, ycen, hix, hiy, ang); + imellipsei(g, 0, xsno+1, flag, type, x, y, xcen, ycen, lox, loy, ang); + for(i=0; i<n; i++){ + imellipsei(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lox+((i+1)*dincx), loy+((i+1)*dincy), ang); + } +} + +void imnpiei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double lo, double hi, int n) +{ + int i; + int xsno; + double dinc; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + while( lo > hi ) lo -= 360.0; + dinc = (hi - lo)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init all shapes */ + impiei(g, 0, xsno, flag, type, x, y, xcen, ycen, lo, hi); + for(i=0; i<n; i++){ + impiei(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lo+(i*dinc), lo+((i+1)*dinc)); + } +} + +void impandai(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double anglo, double anghi, double angn, + double radlo, double radhi, double radn) +{ + int a, r; + int ahi, rhi; + int xsno; + int n=0; + double ainc, rinc; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + while( anglo > anghi ) anglo -= 360.0; + ainc = (anghi - anglo)/angn; + ahi = (int)angn; + rinc = (radhi - radlo)/radn; + rhi = (int)radn; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init pies and annuli */ + imannulusi(g, 0, xsno, flag, type, x, y, xcen, ycen, radlo, radhi); + imqtpiei(g, 0, xsno+1, flag, type, x, y, xcen, ycen, anglo, anghi); + for(a=0; a<ahi; a++){ + for(r=0; r<rhi; r++){ + imannulusi(g, rno+n, sno+(2*n), flag, type, x, y, + xcen, ycen, radlo+(r*rinc), radlo+((r+1)*rinc)); + imqtpiei(g, rno+n, sno+(2*n+1), flag, type, x, y, + xcen, ycen, anglo+(a*ainc), anglo+((a+1)*ainc)); + n++; + } + } +} + +void imbpandai(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double anglo, double anghi, double angn, + double xlo, double ylo, double xhi, double yhi, double radn, + double ang) +{ + int a, r; + int ahi, rhi; + int xsno; + int n=0; + double ainc, xinc, yinc; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + anglo += ang; + anghi += ang; + while( anglo > anghi ) anglo -= 360.0; + ainc = (anghi - anglo)/angn; + ahi = (int)angn; + xinc = (xhi - xlo)/radn; + yinc = (yhi - ylo)/radn; + rhi = (int)radn; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init pies and ellipses */ + imboxi(g, 0, xsno, flag, type, x, y, xcen, ycen, xhi, yhi, ang); + imqtpiei(g, 0, xsno+1, flag, type, x, y, xcen, ycen, anglo, anghi); + imboxi(g, 0, xsno+2, flag, type, x, y, xcen, ycen, xlo, ylo, ang); + for(a=0; a<ahi; a++){ + for(r=1; r<=rhi; r++){ + imboxi(g, rno+n, sno+(2*n), flag, type, x, y, + xcen, ycen, xlo+(r*xinc), ylo+(r*yinc), ang); + imqtpiei(g, rno+n, sno+(2*n+1), flag, type, x, y, + xcen, ycen, anglo+(a*ainc), anglo+((a+1)*ainc)); + n++; + } + } +} + +void imepandai(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double anglo, double anghi, double angn, + double xlo, double ylo, double xhi, double yhi, double radn, + double ang) +{ + int a, r; + int ahi, rhi; + int xsno; + int n=0; + double ainc, xinc, yinc; + + /* NB: do not use x and y variables, they have bogus values */ + + /* get limits */ + anglo += ang; + anghi += ang; + while( anglo > anghi ) anglo -= 360.0; + ainc = (anghi - anglo)/angn; + ahi = (int)angn; + xinc = (xhi - xlo)/radn; + yinc = (yhi - ylo)/radn; + rhi = (int)radn; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + /* init pies and ellipses */ + imellipsei(g, 0, xsno, flag, type, x, y, xcen, ycen, xhi, yhi, ang); + imqtpiei(g, 0, xsno+1, flag, type, x, y, xcen, ycen, anglo, anghi); + imellipsei(g, 0, xsno+2, flag, type, x, y, xcen, ycen, xlo, ylo, ang); + for(a=0; a<ahi; a++){ + for(r=1; r<=rhi; r++){ + imellipsei(g, rno+n, sno+(2*n), flag, type, x, y, + xcen, ycen, xlo+(r*xinc), ylo+(r*yinc), ang); + imqtpiei(g, rno+n, sno+(2*n+1), flag, type, x, y, + xcen, ycen, anglo+(a*ainc), anglo+((a+1)*ainc)); + n++; + } + } +} + +/* accelerator regions -- lower and upper limits are given for n regions */ + +int imnannulus(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, + double lo, double hi, int n) +{ + int i; + int xsno; + double dinc; + + /* get limits */ + dinc = (hi - lo)/(double)n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imannulus(g, 0, xsno, flag, type, x, y, xcen, ycen, lo, hi) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=0; i<n; i++){ + if( imannulus(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lo+(i*dinc), lo+((i+1)*dinc)) ){ + return(1); + } + } + return(0); + } + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imannulus(g, 0, xsno, 1, type, x, y, xcen, ycen, lo, hi) ){ + return(1); + } + return(0); + } +} + +int imnbox(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, + double lox, double loy, double hix, double hiy, int n, + double ang) +{ + int i; + int xsno; + double dincx; + double dincy; + + /* get limits */ + dincx = (hix - lox)/n; + dincy = (hiy - loy)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imbox(g, 0, xsno, flag, type, x, y, xcen, ycen, hix, hiy, ang) ){ + return(0); + } + /* if its in the inner region we lose */ + if( imbox(g, 0, xsno+1, flag, type, x, y, xcen, ycen, lox, loy, ang) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=0; i<n; i++){ + if( imbox(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lox+((i+1)*dincx), loy+((i+1)*dincy), ang) ){ + return(1); + } + } + return(0); + } + /* for excludes, we have to check that we are not in any of them */ + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imbox(g, 0, xsno, 1, type, x, y, xcen, ycen, hix, hiy, ang) ){ + return(1); + } + /* if its in the inner region we win */ + if( imbox(g, 0, xsno+1, 1, type, x, y, xcen, ycen, lox, loy, ang) ){ + return(1); + } + return(0); + } +} + +int imnellipse(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, + double lox, double loy, double hix, double hiy, int n, + double ang) +{ + int i; + int xsno; + double dincx; + double dincy; + + /* get limits */ + dincx = (hix - lox)/n; + dincy = (hiy - loy)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imellipse(g, 0, xsno, flag, type, x, y, xcen, ycen, hix, hiy, ang) ){ + return(0); + } + /* if its in the inner region we lose */ + if( imellipse(g, 0, xsno+1, flag, type, x, y, xcen, ycen, lox, loy, ang) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=0; i<n; i++){ + if( imellipse(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lox+((i+1)*dincx), loy+((i+1)*dincy), ang) ){ + return(1); + } + } + return(0); + } + /* for excludes, we have to check that we are not in any of them */ + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imellipse(g, 0, xsno, 1, type, x, y, xcen, ycen, hix, hiy, ang) ){ + return(1); + } + /* if its in the inner region we win */ + if( imellipse(g, 0, xsno+1, 1, type, x, y, xcen, ycen, lox, loy, ang) ){ + return(1); + } + return(0); + } +} + +int imnpie(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, + double lo, double hi, int n) +{ + int i; + int xsno; + double dinc; + + /* get limits */ + while( lo > hi ) lo -= 360.0; + dinc = (hi - lo)/n; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* check limits */ + if( !impie(g, 0, xsno, flag, type, x, y, xcen, ycen, lo, hi) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=0; i<n; i++){ + if( impie(g, rno+i, sno+i, flag, type, x, y, + xcen, ycen, lo+(i*dinc), lo+((i+1)*dinc)) ){ + return(1); + } + } + return(0); + } + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !impie(g, 0, xsno, 1, type, x, y, xcen, ycen, lo, hi) ){ + return(1); + } + return(0); + } +} + +int impanda(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double anglo, double anghi, double angn, + double radlo, double radhi, double radn) +{ + + int a, r; + int ahi, rhi; + int xsno; + int n=0; + double ainc, rinc; + + /* get limits */ + while( anglo > anghi ) anglo -= 360.0; + ainc = (anghi - anglo)/angn; + ahi = (int)angn; + rinc = (radhi - radlo)/radn; + rhi = (int)radn; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imannulus(g, 0, xsno, flag, type, x, y, xcen, ycen, radlo, radhi) || + !impie(g, 0, xsno+1, flag, type, x, y, xcen, ycen, anglo, anghi) ){ + return(0); + } + /* look through all of them to find the right one */ + for(a=1; a<=ahi; a++){ + for(r=1; r<=rhi; r++){ + if( imannulus(g, rno+n, sno+(2*n), flag, type, x, y, + xcen, ycen, radlo+((r-1)*rinc), radlo+(r*rinc)) && + impie(g, rno+n, sno+(2*n+1), flag, type, x, y, + xcen, ycen, anglo+((a-1)*ainc), anglo+(a*ainc)) ){ + return(1); + } + n++; + } + } + return(0); + } + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imannulus(g, 0, xsno, 1, type, x, y, xcen, ycen, radlo, radhi) ) + return(1); + else if( !impie(g, 0, xsno+1, 1, type, x, y, xcen, ycen, anglo, anghi) ){ + return(1); + } + else{ + return(0); + } + } +} + +int imbpanda(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double anglo, double anghi, double angn, + double xlo, double ylo, double xhi, double yhi, double radn, + double ang) +{ + + int a, r; + int ahi, rhi; + int xsno; + int n=0; + double ainc, xinc, yinc; + + /* get limits */ + anglo += ang; + anghi += ang; + while( anglo > anghi ) anglo -= 360.0; + ainc = (anghi - anglo)/angn; + ahi = (int)angn; + xinc = (xhi - xlo)/radn; + yinc = (yhi - ylo)/radn; + rhi = (int)radn; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imbox(g, 0, xsno, flag, type, x, y, xcen, ycen, xhi, yhi, + ang) ){ + return(0); + } + /* but if its in the inner region we lose */ + else if( imbox(g, 0, xsno+2, flag, type, x, y, xcen, ycen, xlo, ylo, + ang) ){ + return(0); + } + /* its in the box .. must also be in the pie */ + else if( !impie(g, 0, xsno+1, flag, type, x, y, xcen, ycen, anglo, anghi)){ + return(0); + } + /* look through all of them to find the right one */ + for(a=0; a<ahi; a++){ + for(r=1; r<=rhi; r++){ + if( imbox(g, rno+n, sno+(2*n), flag, type, x, y, + xcen, ycen, xlo+(r*xinc), ylo+(r*yinc), ang) && + imqtpie(g, rno+n, sno+(2*n+1), flag, type, x, y, + xcen, ycen, anglo+(a*ainc), anglo+((a+1)*ainc)) ){ + return(1); + } + n++; + } + } + return(0); + } + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imbox(g, 0, xsno, 1, type, x, y, xcen, ycen, xhi, yhi, ang) ) + return(1); + /* if its in the inner region we win */ + else if( !imbox(g, 0, xsno+2, 1, type, x, y, xcen, ycen, xlo, ylo, + ang) ) + return(1); + /* if its not in the pie, we win */ + else if( !impie(g, 0, xsno+1, 1, type, x, y, xcen, ycen, anglo, anghi) ){ + return(1); + } + else{ + return(0); + } + } +} + +int imepanda(GFilt g, int rno, int sno, int flag, int type, + double x, double y, + double xcen, double ycen, + double anglo, double anghi, double angn, + double xlo, double ylo, double xhi, double yhi, double radn, + double ang) +{ + + int a, r; + int ahi, rhi; + int xsno; + int n=0; + double ainc, xinc, yinc; + + /* get limits */ + anglo += ang; + anghi += ang; + while( anglo > anghi ) anglo -= 360.0; + ainc = (anghi - anglo)/angn; + ahi = (int)angn; + xinc = (xhi - xlo)/radn; + yinc = (yhi - ylo)/radn; + rhi = (int)radn; + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imellipse(g, 0, xsno, flag, type, x, y, xcen, ycen, xhi, yhi, + ang) ){ + return(0); + } + /* but if its in the inner region we lose */ + else if( imellipse(g, 0, xsno+2, flag, type, x, y, xcen, ycen, xlo, ylo, + ang) ){ + return(0); + } + /* its in the ellipse .. must also be in the pie */ + else if( !impie(g, 0, xsno+1, flag, type, x, y, xcen, ycen, anglo, anghi)){ + return(0); + } + /* look through all of them to find the right one */ + for(a=0; a<ahi; a++){ + for(r=1; r<=rhi; r++){ + if( imellipse(g, rno+n, sno+(2*n), flag, type, x, y, + xcen, ycen, xlo+(r*xinc), ylo+(r*yinc), ang) && + imqtpie(g, rno+n, sno+(2*n+1), flag, type, x, y, + xcen, ycen, anglo+(a*ainc), anglo+((a+1)*ainc)) ){ + return(1); + } + n++; + } + } + return(0); + } + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imellipse(g, 0, xsno, 1, type, x, y, xcen, ycen, xhi, yhi, ang) ) + return(1); + /* if its in the inner region we win */ + else if( !imellipse(g, 0, xsno+2, 1, type, x, y, xcen, ycen, xlo, ylo, + ang) ) + return(1); + /* if its not in the pie, we win */ + else if( !impie(g, 0, xsno+1, 1, type, x, y, xcen, ycen, anglo, anghi) ){ + return(1); + } + else{ + return(0); + } + } +} + +#ifdef __STDC__ +void +imvannulusi(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, n; + int maxpts; + int xsno; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvannulusi(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double ang; + double *xv; + int i, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + /* this should be impossible ... */ + if( n == 2 ){ + imannulusi(g, rno, sno, flag, type, x, y, xcen, ycen, xv[0], xv[1]); + return; + } + imannulusi(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[0], xv[n-1]); + for(i=0; i<(n-1); i++){ + imannulusi(g, rno+i, sno+i, flag, type, x, y, xcen, ycen, xv[i], xv[i+1]); + } +} + +#ifdef __STDC__ +void +imvboxi(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, j, n; + int maxpts; + int xsno; + double ang; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvboxi(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double ang; + double *xv; + int i, j, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + ang = xv[--n]; + /* this should be impossible ... */ + if( n == 2 ){ + imboxi(g, rno, sno, flag, type, x, y, xcen, ycen, xv[0], xv[1], ang); + return; + } + imboxi(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[n-2], xv[n-1], ang); + imboxi(g, 0, xsno+1, flag, type, x, y, xcen, ycen, xv[0], xv[1], ang); + for(i=2, j=0; i<n; i+=2, j++){ + imboxi(g, rno+j, sno+j, flag, type, x, y, xcen, ycen, xv[i], xv[i+1], ang); + } +} + +#ifdef __STDC__ +void +imvellipsei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, j, n; + int maxpts; + int xsno; + double ang; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvellipsei(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double ang; + double *xv; + int i, j, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + ang = xv[--n]; + /* this should be impossible ... */ + if( n == 2 ){ + imellipsei(g, rno, sno, flag, type, x, y, xcen, ycen, xv[0], xv[1], ang); + return; + } + imellipsei(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[n-2], xv[n-1], ang); + imellipsei(g, 0, xsno+1, flag, type, x, y, xcen, ycen, xv[0], xv[1], ang); + for(i=2, j=0; i<n; i+=2, j++){ + imellipsei(g, rno+j, sno+j, flag, type, x, y, xcen, ycen, + xv[i], xv[i+1], ang); + } +} + +#ifdef __STDC__ +void +imvpiei(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, n; + int maxpts; + int xsno; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvpiei(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double ang; + double *xv; + int i, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + /* this should be impossible ... */ + if( n == 2 ){ + impiei(g, rno, sno, flag, type, x, y, xcen, ycen, xv[0], xv[1]); + return; + } + impiei(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[0], xv[n-1]); + for(i=0; i<(n-1); i++){ + impiei(g, rno+i, sno+i, flag, type, x, y, xcen, ycen, xv[i], xv[i+1]); + } +} + +#ifdef __STDC__ +void +imvpointi(GFilt g, int rno, int sno, int flag, int type, + double x, double y, ...) +{ + int i, j, n; + int maxpts; + int xsno; + double *xv; + va_list args; + va_start(args, y); +#else +int imvpointi(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double *xv; + int i, j, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + for(i=0, j=0; i<n; i+=2, j++){ + impointi(g, rno+j, sno+j, flag, type, x, y, xv[i], xv[i+1]); + } +} + +/* varargs regions -- a series of lower and upper limits is specified */ + +#ifdef __STDC__ +int +imvannulus(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, n; + int maxpts; + int xsno; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvannulus(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double *xv; + int i, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + /* this should be impossible ... */ + if( n == 2 ){ + return(imannulus(g, rno, sno, flag, type, x, y, xcen, ycen, xv[0], xv[1])); + } + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imannulus(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[0], xv[n-1]) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=0; i<n; i++){ + if( imannulus(g, rno+i, sno+i, flag, type, x, y, xcen, ycen, + xv[i], xv[i+1]) ){ + return(1); + } + } + return(0); + } + /* for excludes, we have to check that we are not in any of them */ + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imannulus(g, 0, xsno, 1, type, x, y, xcen, ycen, xv[0], xv[n-1]) ){ + return(1); + } + return(0); + } +} + +#ifdef __STDC__ +int +imvbox(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, j, n; + int maxpts; + int xsno; + double ang; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvbox(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double ang; + double *xv; + int i, j, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + ang = xv[--n]; + /* this should be impossible ... */ + if( n == 2 ){ + return(imbox(g, rno, sno, flag, type, x, y, + xcen, ycen, xv[0], xv[1], ang)); + } + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imbox(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[n-2], xv[n-1], + ang) ){ + return(0); + } + /* if its in the inner region we lose */ + if( imbox(g, 0, xsno+1, flag, type, x, y, xcen, ycen, xv[0], xv[1], ang) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=2, j=0; i<n; i+=2, j++){ + if( imbox(g, rno+j, sno+j, flag, type, x, y, xcen, ycen, + xv[i], xv[i+1], ang) ){ + return(1); + } + } + return(0); + } + /* for excludes, we have to check that we are not in any of them */ + else{ + /* if its not somewhere inside the entire region we win ... */ + if( !imbox(g, 0, xsno, 1, type, x, y, xcen, ycen, xv[n-2], xv[n-1], ang) ){ + return(1); + } + /* if its in the inner region we win */ + else if( imbox(g, 0, xsno+1, 1, type, x, y, xcen, ycen, xv[0], xv[1], ang) ){ + return(1); + } + return(0); + } +} + +#ifdef __STDC__ +int +imvellipse(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, j, n; + int maxpts; + int xsno; + double ang; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvellipse(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double ang; + double *xv; + int i, j, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + ang = xv[--n]; + /* this should be impossible ... */ + if( n == 2 ){ + return(imellipse(g, rno, sno, flag, type, x, y, + xcen, ycen, xv[0], xv[1], ang)); + } + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !imellipse(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[n-2], xv[n-1], + ang) ){ + return(0); + } + /* if its in the inner region we lose */ + if( imellipse(g, 0, xsno+1, flag, type, x, y, xcen, ycen, xv[0], xv[1], + ang) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=2, j=0; i<n; i+=2, j++){ + if( imellipse(g, rno+j, sno+j, flag, type, x, y, xcen, ycen, + xv[i], xv[i+1], ang) ){ + return(1); + } + } + return(0); + } + /* for excludes, we have to check that we are not in any of them */ + else{ + /* if its not somewhere inside the entire region we lose ... */ + if( !imellipse(g, 0, xsno, 1, type, x, y, xcen, ycen, xv[n-2], xv[n-1], + ang) ){ + return(1); + } + /* if its in the inner region we win */ + if( imellipse(g, 0, xsno+1, 1, type, x, y, xcen, ycen, xv[0], xv[1], ang) ){ + return(1); + } + return(0); + } +} + +#ifdef __STDC__ +int +imvpie(GFilt g, int rno, int sno, int flag, int type, + double x, double y, double xcen, double ycen, ...) +{ + int i, n; + int maxpts; + int xsno; + double *xv; + va_list args; + va_start(args, ycen); +#else +int imvpie(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double xcen, ycen; + double *xv; + int i, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); + xcen = va_arg(args, double); + ycen = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + /* this should be impossible ... */ + if( n == 2 ){ + return(impie(g, rno, sno, flag, type, x, y, + xcen, ycen, xv[0], xv[1])); + } + if( flag ){ + /* if its not somewhere inside the entire region we lose ... */ + if( !impie(g, 0, xsno, flag, type, x, y, xcen, ycen, xv[0], xv[n-1]) ){ + return(0); + } + /* look through all of them to find the right one */ + for(i=0; i<n; i++){ + if( impie(g, rno+i, sno+i, flag, type, x, y, xcen, ycen, xv[i], xv[i+1]) ){ + return(1); + } + } + return(0); + } + /* for excludes, we have to check that we are not in any of them */ + else{ + /* if its not somewhere inside the entire region we lose ... */ + if( !impie(g, 0, xsno, 1, type, x, y, xcen, ycen, xv[0], xv[n-1]) ){ + return(1); + } + return(1); + } +} + +#ifdef __STDC__ +int +imvpoint(GFilt g, int rno, int sno, int flag, int type, + double x, double y, ...) +{ + int i, j, n; + int maxpts; + int xsno; + double *xv; + va_list args; + va_start(args, y); +#else +int imvpoint(va_alist) va_dcl +{ + GFilt g; + int rno, sno, flag, type; + double x, y; + double *xv; + int i, j, n; + int maxpts; + int xsno; + va_list args; + va_start(args); + g = va_arg(args, GFilt); + rno = va_arg(args, int); + sno = va_arg(args, int); + flag = va_arg(args, int); + type = va_arg(args, int); + x = va_arg(args, double); + y = va_arg(args, double); +#endif + xsno = (g->nshapes+1)+((sno-1)*XSNO); + if( !g->shapes[xsno].xv ){ + maxpts = MASKINC; + g->shapes[xsno].xv = (double *)calloc(maxpts, sizeof(double)); + g->shapes[xsno].nv = 0; + while( 1 ){ + if( g->shapes[xsno].nv >= maxpts ){ + maxpts += MASKINC; + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + maxpts*sizeof(double)); + } + g->shapes[xsno].xv[g->shapes[xsno].nv] = va_arg(args, double); + if( feq(g->shapes[xsno].xv[g->shapes[xsno].nv],PSTOP) && + feq(g->shapes[xsno].xv[g->shapes[xsno].nv-1],PSTOP) ){ + g->shapes[xsno].nv--; + break; + } + g->shapes[xsno].nv++; + } + va_end(args); + g->shapes[xsno].xv = (double *)realloc(g->shapes[xsno].xv, + g->shapes[xsno].nv*sizeof(double)); + } + n = g->shapes[xsno].nv; + xv = g->shapes[xsno].xv; + /* look through all of them to find the right one */ + for(i=0, j=0; i<n; i+=2, j++){ + if( impoint(g, rno+j, sno+j, flag, type, x, y, xv[i], xv[i+1]) ){ + return(1); + } + } + return(0); +} + +void imimagemaski(GFilt g, int UNUSED(rno), int UNUSED(sno), + int UNUSED(flag), int UNUSED(type), + double UNUSED(x), double UNUSED(y)) +{ + int i, j; + int got; + int fsize; + int mblock; + FilterMask xmasks; + + /* NB: do not use x and y variables, they have bogus values */ + /* scale mask x, y values to match the incoming image section */ + mblock = (int)(((double)((g->xmax - g->xmin + 1)/g->block)/(double)g->maskdim) + 0.5); + if( mblock < 1 ){ + g->masks = NULL; + g->nmask = 0; + } + /* copy all masks and translate x,y positions as needed */ + fsize = g->nmask * sizeof(FilterMaskRec) * mblock; + xmasks = (FilterMask)calloc(1, fsize); + for(got=0, i=0; i<g->nmask; i++){ + xmasks[got].region = g->masks[i].region; + xmasks[got].y = (g->masks[i].y - 1.0) * mblock + 1.0; + xmasks[got].xstart = (g->masks[i].xstart - 1.0) * mblock + 1.0; + xmasks[got].xstop = (g->masks[i].xstop - 1.0) * mblock + 1.0; + /* replicate the segment up to the block factor */ + for(j=1; j<mblock; j++){ + xmasks[got+j].region = xmasks[got].region; + xmasks[got+j].y = xmasks[got].y+j; + xmasks[got+j].xstart = xmasks[got].xstart; + xmasks[got+j].xstop = xmasks[got].xstop; + } + got += mblock; + } + /* sort by y and x */ + qsort(xmasks, got, sizeof(FilterMaskRec), imagemaskcmp); + /* we don't free previous masks, just overwrite, because they were not + allocated in the first place */ + /* if( g->masks ) free(g->masks); */ + /* set up new mask records */ + g->masks = xmasks; + g->nmask = got; + + /* now mark each y line that has a y mask value */ + for(i=0; i<got; i++){ + if( !g->ybuf[g->masks[i].y] ) + g->ybuf[g->masks[i].y] = i+1; + } +} + +int imimagemask(GFilt g, int UNUSED(rno), int UNUSED(sno), + int UNUSED(flag), int UNUSED(type), + double x, double y) +{ + int i; + int ix, iy; + + if( g->nmask == 0 ) + return(0); + if( g->evsect ){ + if( g->usebinsiz ){ + x = (int)((x - g->tlminx)/g->binsizx + 1.0); + y = (int)((y - g->tlminy)/g->binsizy + 1.0); + } + else{ + x = (int)((x - g->tlminx) + 1.0); + y = (int)((y - g->tlminy) + 1.0); + } + x = (int)((x - g->xmin)/g->block + 1.0); + y = (int)((y - g->ymin)/g->block + 1.0); + if(y < g->y0) return 0; + if(y > g->y1) return 0; + if(x < g->x0) return 0; + if(x > g->x1) return 0; + } + ix = (int)x; + iy = (int)y; + /* look for mask segment containing ix, iy */ + i = g->ybuf[iy]; + /* ybuf contains offset+1 of first mask record containing iy */ + if( i != 0 ){ + i--; + for(; i<g->nmask; i++){ + if( g->masks[i].y > iy ){ + return(0); + } + if( (ix >= g->masks[i].xstart) && (ix <= g->masks[i].xstop) ){ + g->rid = g->masks[i].region; + return(1); + } + } + } + return(0); +} + |