summaryrefslogtreecommitdiffstats
path: root/funtools/funcen.c
diff options
context:
space:
mode:
Diffstat (limited to 'funtools/funcen.c')
-rw-r--r--funtools/funcen.c830
1 files changed, 830 insertions, 0 deletions
diff --git a/funtools/funcen.c b/funtools/funcen.c
new file mode 100644
index 0000000..a2efae0
--- /dev/null
+++ b/funtools/funcen.c
@@ -0,0 +1,830 @@
+/*
+ * Copyright (c) 2004 Smithsonian Astrophysical Observatory
+ */
+
+/*
+ * osx, linux:
+ * gcc -O2 -g -o funcen funcen.c -I $HOME/include -L$HOME/lib -lfuntools -lm
+ *
+ * solaris:
+ * gcc -O2 -g -o funcen funcen.c -I /soft/saord/include -L/soft/saord/lib \
+ * -lfuntools -lsocket -lnsl -ldl -lm
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <funtools.h>
+#include <filter.h>
+#include <swap.h>
+#include <word.h>
+#include <wcs.h>
+#include <xalloc.h>
+
+#define MAXROW 8192
+#define MAXREGIONS 8192
+#define WCSSYS "j2000"
+
+extern char *optarg;
+extern int optind;
+
+/* single event (used in array of structs) */
+typedef struct evstruct{
+ int region;
+ double x, y;
+} *Ev, EvRec;
+
+#ifdef ANSI_FUNC
+static void
+usage (char *fname)
+#else
+static void usage(fname)
+ char *fname;
+#endif
+{
+ fprintf(stderr,
+ "usage: %s <switches> eventfile region(s)\n", fname);
+ fprintf(stderr, "\n");
+ fprintf(stderr, "optional switches:\n");
+ fprintf(stderr, " -i\t\t# use image filtering (default: event filtering)\n");
+ fprintf(stderr, " -n iter\t# max number of iterations (default: 0)\n");
+ fprintf(stderr, " -t tol\t# pixel tolerance distance (default: 1.0)\n");
+ fprintf(stderr, " -v [0,1,2,3]\t# output verbosity level (default: 0)\n");
+ fprintf(stderr, "\n");
+ fprintf(stderr,
+ "The default output consists of one line for each region containing:\n");
+ fprintf(stderr, "cnts x y [ra dec]. Use -v for prettier output and diagnostics.\n");
+ fprintf(stderr,
+ "Filter extensions are permitted as part of the file specification.\n");
+ fprintf(stderr,
+ "WCS info, if available, matches the coordinate system of the data file.\n");
+ fprintf(stderr,
+ "Input regions can be in physical or wcs coordinates (with coordsys),\n");
+ fprintf(stderr,
+ "Multiple regions and region files are supported but compound regions,\n");
+ fprintf(stderr, "polygons, and points are not permitted.\n");
+ fprintf(stderr, "\n");
+ fprintf(stderr, "examples:\n");
+ fprintf(stderr,
+ " %s -n 5 -t .1 snr.ev'[pha=1:3]' 'cir(512,512,10);cir 450 450 10'\n",
+ fname);
+ fprintf(stderr,
+ " %s -n 5 -t .1 snr.ev 'fk5;circle(345.292335 58.881625,78.458952\")'\n",
+ fname);
+ fprintf(stderr, "\n(version: %s)\n", FUN_VERSION);
+ exit(1);
+}
+
+#ifdef ANSI_FUNC
+int
+main (int argc, char **argv)
+#else
+int
+main(argc, argv)
+ int argc;
+ char **argv;
+#endif
+{
+ int c;
+ int iter;
+ int i;
+ int got;
+ int ntol;
+ int len;
+ int dowcs=0;
+ int doimfilt=0;
+ int domem=0;
+ int ip=0, wcnt=0, incl=0;
+ int verbose=0;
+ int args=0;
+ int maxrow=MAXROW;
+ int niter=0;
+ int nregion=0;
+ int type, endian, x0, x1, y0, y1, block, rawsize;
+ int nev[MAXREGIONS+1];
+ int tflag[MAXREGIONS+1];
+ int *rbuf=NULL;
+ int cfd;
+ int clen=0;
+ int tltyp[2];
+ int bin[2];
+ double tlmin[2];
+ // double tlmax[2];
+ double binsiz[2];
+ double xtot[MAXREGIONS+1], ytot[MAXREGIONS+1];
+ double xavg[MAXREGIONS+1], yavg[MAXREGIONS+1];
+ double ox[MAXREGIONS+1], oy[MAXREGIONS+1];
+ double pxavg[MAXREGIONS+1], pyavg[MAXREGIONS+1];
+ double xim, yim;
+ double ra, dec;
+ double ctolsq=0.0, tolsq=0.0;
+ double tol=1.0;
+ char *s;
+ char *fname=NULL;
+ char *iregion=NULL;
+ char *cname=NULL;
+ char *cregion=NULL;
+ char *rawbuf=NULL;
+ char *bincols=NULL;
+ char *tmpdir=NULL;
+ char prefix[SZ_LINE];
+ char cregionfile[SZ_LINE];
+ char mode[SZ_LINE];
+ char lastd[2];
+ char extn[3];
+ char fname2[SZ_LINE];
+ char coorsys[SZ_LINE];
+ char wcssys[SZ_LINE];
+ char tbuf[SZ_LINE];
+ char tbuf2[SZ_LINE];
+ char *regions[MAXREGIONS+1];
+ struct stat statbuf;
+ FILE *fd;
+ Fun fun=NULL;
+ Ev ev, ebuf=NULL;
+ FITSHead header=NULL;
+ Filter filter=NULL;
+ struct WorldCoor *wcs=NULL;
+
+ /* exit on gio errors */
+ if( !getenv("GERROR") )
+ setgerror(2);
+
+ /* we want the args in the same order in which they arrived, and
+ gnu getopt sometimes changes things without this */
+ putenv("POSIXLY_CORRECT=true");
+
+ /* init wcssys to default */
+ strncpy(wcssys, WCSSYS, SZ_LINE-1);
+
+ /* process switch arguments */
+ while ((c = getopt(argc, argv, "imn:t:v:w:")) != -1){
+ switch(c){
+ case 'i':
+ doimfilt=1;
+ break;
+ case 'm':
+ domem=1;
+ break;
+ case 'n':
+ niter = atoi(optarg);
+ if( niter < 0 ) niter = 0;
+ break;
+ case 't':
+ tol = strtod(optarg, (char **)NULL);
+ break;
+ case 'v':
+ verbose = atoi(optarg);
+ if( verbose < 0 ) verbose = 0;
+ break;
+ case 'w':
+ strncpy(wcssys, optarg, SZ_LINE-1);
+ break;
+ }
+ }
+
+ /* get maxrow, if user-specified */
+ if( (s=(char *)getenv("FUN_MAXROW")) != NULL )
+ maxrow = atoi(s);
+
+ /* check for required arguments */
+ args = argc - optind;
+ if( args < 2 ) usage(argv[0]);
+ /* arg 1: input file name */
+ fname = argv[optind+0];
+ /* arg 2: initial region */
+ iregion = argv[optind+1];
+
+ /* intialization: clear buffers */
+ memset(regions, 0, sizeof(char *)*(MAXREGIONS+1));
+ memset(ox, 0, sizeof(double)*(MAXREGIONS+1));
+ memset(oy, 0, sizeof(double)*(MAXREGIONS+1));
+ memset(pxavg, 0, sizeof(double)*(MAXREGIONS+1));
+ memset(pyavg, 0, sizeof(double)*(MAXREGIONS+1));
+ memset(tflag, 0, sizeof(int)*(MAXREGIONS+1));
+ /* deal with the square of the tolerance */
+ tolsq = tol*tol;
+ /* no coords system yet */
+ *coorsys = '\0';
+ /* no temp region file yet */
+ *cregionfile = '\0';
+
+ /* get prefix for filter region temp file */
+ if( !(tmpdir = (char *)getenv("FILTER_TMPDIR")) &&
+ !(tmpdir = (char *)getenv("TMPDIR")) &&
+ !(tmpdir = (char *)getenv("TMP")) )
+ tmpdir = DEFAULT_FILTER_TMPDIR;
+ if( !*tmpdir )
+ tmpdir = ".";
+ snprintf(prefix, SZ_LINE-1, "%s/funcen", tmpdir);
+
+ /* if processing in memory, read all events now */
+ if( domem ){
+ /* open the data file once */
+ if( !(fun = FunOpen(fname, "r", NULL)) ){
+ gerror(stderr, "can't FunOpen file (or find extension): %s\n", fname);
+ }
+ /* make sure we have a table */
+ FunInfoGet(fun, FUN_TYPE, &type, 0);
+ if( type != FUN_TABLE ){
+ gerror(stderr,
+ "please specify a table (image support not yet implemented)\n");
+ }
+ /* for wcs output, get info to allow conversion of phyiscal to image */
+ FunInfoGet(fun, FUN_BINOFFS, bin, 0);
+ for(i=0; i<=1; i++){
+ if( bin[i] >=0 ){
+ FunColumnLookup(fun, NULL, bin[i], NULL,
+ &(tltyp[i]), NULL, NULL, NULL, NULL);
+ tlmin[i] = FunParamGetd(fun, "TLMIN", bin[i]+1, 0.0, &got);
+ // tlmax[i] = FunParamGetd(fun, "TLMAX", bin[i]+1, 0.0, &got);
+ binsiz[i] = FunParamGetd(fun, "TDBIN", bin[i]+1, 1.0, &got);
+ }
+ }
+ /* get row number that we need to read */
+ maxrow = FunParamGeti(fun, "NAXIS", 2, -1, &got);
+ if( (maxrow < 0) || !got ){
+ gerror(stderr, "can't determine total number of events in table\n");
+ }
+ /* select binning columns */
+ got = FunColumnSelect(fun, sizeof(EvRec), NULL,
+ "$x", "D", "r", FUN_OFFSET(Ev, x),
+ "$y", "D", "r", FUN_OFFSET(Ev, y),
+ "$region", "J", "r", FUN_OFFSET(Ev, region),
+ NULL);
+ if( !(ebuf=(Ev)xmalloc(maxrow*sizeof(EvRec))) ){
+ gerror(stderr, "can't allocate memory for all events in table\n");
+ }
+ ebuf = (Ev)FunTableRowGet(fun, ebuf, maxrow, NULL, &got);
+ if( got != maxrow ){
+ gerror(stderr,
+ "could not read all rows into memory (wanted %d, got %d)\n",
+ maxrow, got);
+ }
+ }
+
+ /* if we are filtering internally, allocate the return buffer */
+ if( domem || doimfilt ){
+ if( !(rbuf = xmalloc(maxrow*sizeof(int))) ){
+ gerror(stderr, "can't allocate memory\n");
+ }
+ }
+ /* not filtering internally:, prepare the ifile for new filter each time */
+ else{
+ /* make up the beginning of the file name -- we will append the region */
+ strncpy(fname2, fname, SZ_LINE-2);
+ i = strlen(fname2)-1;
+ /* if we have an extension, prepare to append to it */
+ if( fname2[i] == ']' ){
+ fname2[i] = '\0';
+ strncat(fname2, ",(", SZ_LINE-1);
+ strcpy(extn, ")]");
+ }
+ /* otherwise, start a new extension */
+ else{
+ strncat(fname2, "[", SZ_LINE-1);
+ strcpy(extn, "]");
+ }
+ }
+
+ /* check for region file and retrieve contents, if necessary */
+ if( *iregion != '@' ){
+ clen = SZ_LINE;
+ if( !(cregion = xcalloc(clen, sizeof(char))) ){
+ gerror(stderr, "can't allocate memory\n");
+ }
+ if( !(cname = xcalloc(clen, sizeof(char))) ){
+ gerror(stderr, "can't allocate memory\n");
+ }
+ strncpy(cregion, iregion, SZ_LINE-1);
+ }
+ /* yes, its a file */
+ else{
+ if( (fd=fopen(iregion+1, "r")) == NULL ){
+ gerror(stderr, "can't open file: %s\n", iregion);
+ }
+ if( fstat(fileno(fd), &statbuf) <0 ){
+ gerror(stderr, "can't get file information: %s\n", iregion);
+ }
+ clen = statbuf.st_size + SZ_LINE;
+ if( !(cregion=(char *)xcalloc(clen, sizeof(char))) ){
+ gerror(stderr, "can't allocate memory for region file: %s\n", iregion);
+ }
+ if( !(cname=(char *)xcalloc(clen, sizeof(char))) ){
+ gerror(stderr, "can't allocate memory for region file: %s\n", iregion);
+ }
+ while( fgets(tbuf, SZ_LINE, fd) ){
+ if( *tbuf == '#' || *tbuf == '\n' ) continue;
+ if( (s=strchr(tbuf, '#')) ) *s = '\0';
+ nowhite(tbuf, tbuf2);
+ len = strlen(tbuf2);
+ if( len > 1 ){
+ strncat(cregion, tbuf2, SZ_LINE-1);
+ if( tbuf2[len-1] != ';' )
+ strncat(cregion, ";", SZ_LINE-1);
+ }
+ }
+ fclose(fd);
+ }
+
+ /* process each region to make up a region formatting string */
+ newdtable(";\n,()");
+ while( word(cregion, tbuf, &ip) ){
+ lastd[0] = lastdelim();
+ lastd[1] = '\0';
+ switch(++wcnt){
+ case 1:
+ /* handle coorsys specially */
+ if( strstr("j2000|fk5|icrs|b1950|fk4|galactic|ecliptic|linear|image|physical|amplifier|detector", tbuf) ){
+ strncat(coorsys, tbuf, SZ_LINE-1);
+ strncat(coorsys, lastd, SZ_LINE-1);
+ wcnt--;
+ continue;
+ }
+ if( strstr(tbuf, "pol") ){
+ gerror(stderr, "polygons are not allowed\n");
+ }
+ if( strstr(tbuf, "poi") ){
+ gerror(stderr, "points are not allowed\n");
+ }
+ if( *tbuf == '-' ){
+ incl = 0;
+ }
+ else{
+ incl = 1;
+ nregion++;
+ }
+ if( !regions[nregion] ){
+ if( !(regions[nregion] = xcalloc(SZ_LINE, sizeof(char))) ){
+ gerror(stderr, "can't allocate memory\n");
+ }
+ }
+ if( *coorsys ){
+ strncat(regions[nregion], coorsys, SZ_LINE-1);
+ }
+ strncat(regions[nregion], tbuf, SZ_LINE-1);
+ break;
+ case 2:
+ if( incl ){
+ ox[nregion] = strtod(tbuf, (char **)NULL);
+ if( *coorsys && strncmp(coorsys, "physical", 8) )
+ strncat(regions[nregion], "%.2fp", SZ_LINE-1);
+ else
+ strncat(regions[nregion], "%.2f", SZ_LINE-1);
+ }
+ else{
+ strncat(regions[nregion], tbuf, SZ_LINE-1);
+ }
+ break;
+ case 3:
+ if( incl ){
+ oy[nregion] = strtod(tbuf, (char **)NULL);
+ if( *coorsys && strncmp(coorsys, "physical", 8) )
+ strncat(regions[nregion], "%.2fp", SZ_LINE-1);
+ else
+ strncat(regions[nregion], "%.2f", SZ_LINE-1);
+ }
+ else{
+ strncat(regions[nregion], tbuf, SZ_LINE-1);
+ }
+ break;
+ default:
+ strncat(regions[nregion], tbuf, SZ_LINE-1);
+ break;
+ }
+ if( (lastd[0] != ';') && (lastd[0] != '\n') ){
+ strncat(regions[nregion], lastd, SZ_LINE-1);
+ }
+ if( (lastd[0] == ';') || (lastd[0] == ')') || (lastd[0] == '\n') ){
+ wcnt = 0;
+ *coorsys = '\0';
+ }
+ }
+ freedtable();
+
+ /* iterate n times (or until tolerance is reached) to calculate centroid */
+ for(iter=0; iter<niter; iter++){
+ /* reset temp buffers */
+ for(i=1; i<=nregion; i++){
+ if( !tflag[i] ){
+ xtot[i] = 0;
+ ytot[i] = 0;
+ xavg[i] = 0;
+ yavg[i] = 0;
+ nev[i] = 0;
+ }
+ }
+ /* process file: re-open the input file, possibly with a new region */
+ if( !domem ){
+ /* append current region to filename root */
+ if( !doimfilt ){
+ /* if the region string is not long, we append it to the file name */
+ if( (strlen(fname2) + strlen(cregion) + strlen(extn)) < (SZ_LINE-1) ){
+ snprintf(cname, clen-1, "%s%s%s", fname2, cregion, extn);
+ }
+ /* otherwise we have to use a temp file */
+ else{
+ /* delete previous temp region file */
+ if( *cregionfile ) unlink(cregionfile);
+ /* get new temp file name */
+ if( !(cfd=mkrtemp(prefix, ".reg", cregionfile, SZ_LINE, 1)) ){
+ gerror(stderr, "can't open temp region file: %s\n", cregionfile);
+ }
+ else{
+ if( !write(cfd, cregion, strlen(cregion)) ){
+ gerror(stderr, "can't write temp file: %s\n", cregionfile);
+ }
+ close(cfd);
+ }
+ snprintf(cname, clen-1, "%s@%s%s", fname2, cregionfile, extn);
+ }
+ }
+ /* but not for image-style filtering -- we filter the region manually */
+ else{
+ strncpy(cname, fname, clen-1);
+ }
+ if( verbose >= 2 ){
+ if( iter > 0 ) fprintf(stdout, "\n");
+ fprintf(stdout, "iter%d: %s\n", iter, cname);
+ }
+ /* make sure total data file name is less than SZ_LINE */
+ if( strlen(cname) >= SZ_LINE ){
+ gerror(stderr,
+ "region length (%d) > max allowable (4096 - filename_length)\n",
+ strlen(cregion));
+ }
+ /* open the data file with the new region as part of the filter */
+ if( !(fun = FunOpen(cname, "r", NULL)) ){
+ gerror(stderr, "can't FunOpen file (or find extension): %s\n", fname);
+ }
+ /* make sure we have a table */
+ FunInfoGet(fun, FUN_TYPE, &type, 0);
+ if( type != FUN_TABLE ){
+ gerror(stderr,
+ "please specify a table (image support not yet implemented)\n");
+ }
+ /* for wcs output, get info to allow conversion of phyiscal to image */
+ FunInfoGet(fun, FUN_BINOFFS, bin, 0);
+ for(i=0; i<=1; i++){
+ if( bin[i] >=0 ){
+ FunColumnLookup(fun, NULL, bin[i], NULL,
+ &(tltyp[i]), NULL, NULL, NULL, NULL);
+ tlmin[i] = FunParamGetd(fun, "TLMIN", bin[i]+1, 0.0, &got);
+ // tlmax[i] = FunParamGetd(fun, "TLMAX", bin[i]+1, 0.0, &got);
+ binsiz[i] = FunParamGetd(fun, "TDBIN", bin[i]+1, 1.0, &got);
+ }
+ }
+ got = FunColumnSelect(fun, sizeof(EvRec), NULL,
+ "$x", "D", "r", FUN_OFFSET(Ev, x),
+ "$y", "D", "r", FUN_OFFSET(Ev, y),
+ "$region", "J", "r", FUN_OFFSET(Ev, region),
+ NULL);
+ }
+ else{
+ if( verbose >= 2 ){
+ if( iter > 0 ) fprintf(stdout, "\n");
+ fprintf(stdout, "iter%d: %s\n", iter, cregion);
+ }
+ }
+ /* for internal filtering, open up a new filter */
+ if( domem || doimfilt ){
+ FunInfoGet(fun, FUN_ENDIAN, &endian, FUN_HEADER, &header,
+ FUN_SECT_X0, &x0, FUN_SECT_X1, &x1,
+ FUN_SECT_Y0, &y0, FUN_SECT_Y1, &y1,
+ FUN_SECT_BLOCK, &block, FUN_BINCOLS, &bincols,
+ FUN_RAWSIZE, &rawsize,
+ 0);
+ if( doimfilt ){
+ /* make sure we will have a valid image section in which to filter */
+ if( (x1-x0 <= 0) || (y1-y0 <= 0) )
+ gerror(stderr, "invalid or zero image dimensions(s) for table\n",
+ bincols ? bincols : "invalid bincols?");
+ /* make up the mode string */
+ snprintf(mode, SZ_LINE-1,
+ "type=events,convert=%s,evsect=\"%d %d %d %d %d\"",
+ (endian == is_bigendian()) ? "false" : "true",
+ x0, x1, y0, y1, block);
+ }
+ else{
+ snprintf(mode, SZ_LINE-1,
+ "type=events,convert=%s",
+ (endian == is_bigendian()) ? "false" : "true");
+ }
+ /* add columns */
+ if( bincols ){
+ strcat(mode, ",");
+ strcat(mode, bincols);
+ }
+ filter = FilterOpen(header, cregion, mode);
+ }
+
+ /* filter events and accumulate total events for each region */
+ while( 1 ){
+ /* get data (already have it for memory processing */
+ if( !domem ){
+ if( !(ebuf = (Ev)FunTableRowGet(fun, NULL, maxrow, NULL, &got)) )
+ break;
+ }
+ /* do standard event filtering */
+ if( !domem && !doimfilt ){
+ for(i=0; i<got; i++){
+ ev = ebuf+i;
+ if( ev->region > nregion ){
+ gerror(stderr, "compound regions are not allowed!\n");
+ }
+ if( !tflag[ev->region] ){
+ xtot[ev->region] += ev->x;
+ ytot[ev->region] += ev->y;
+ nev[ev->region] += 1;
+ /* debugging output */
+ if( verbose > 2 ){
+ fprintf(stdout, "\t%.2f %.2f %d\n", ev->x, ev->y, ev->region);
+ }
+ }
+ }
+ }
+ /* memory and image-style filtering is done manually */
+ else{
+ /* get pointer to raw buffer */
+ FunInfoGet(fun, FUN_RAWBUF, &rawbuf, 0);
+ /* get events which pass the region filter */
+ if( filter && (filter != NOFILTER) &&
+ FilterEvents(filter, rawbuf, rawsize, got, rbuf)){
+ /* loop through events, process those which are in a region */
+ for(i=0; i<got; i++){
+ if( rbuf[i] > 0 ){
+ ev = ebuf + i;
+ if( !tflag[rbuf[i]] ){
+ xtot[rbuf[i]] += ev->x;
+ ytot[rbuf[i]] += ev->y;
+ nev[rbuf[i]] += 1;
+ }
+ }
+ }
+ }
+ }
+ if( domem ) {
+ break;
+ }
+ else{
+ if( ebuf ) xfree(ebuf);
+ }
+ }
+ /* calculate centroid position for each region */
+ ntol = 0;
+ for(i=1; i<=nregion; i++){
+ if( !tflag[i] ){
+ if( nev[i] > 0 ){
+ xavg[i] = xtot[i] / nev[i];
+ yavg[i] = ytot[i] / nev[i];
+ }
+ else{
+ xavg[i] = 0.0;
+ yavg[i] = 0.0;
+ }
+ /* check tolerance of this centroid */
+ ctolsq = (((xavg[i] - pxavg[i]) * (xavg[i] - pxavg[i])) +
+ ((yavg[i] - pyavg[i]) * (yavg[i] - pyavg[i])));
+ if( ctolsq <= tolsq ){
+ tflag[i] = 1;
+ }
+ else{
+ ntol++;
+ }
+ }
+ if( verbose > 1 ){
+ if( i != 1 ) fprintf(stdout, "\n");
+ snprintf(tbuf, SZ_LINE-1, regions[i], xavg[i], yavg[i]);
+ fprintf(stdout, "\tevents:\t\t%d\n", nev[i]);
+ fprintf(stdout, "\tcentroid(x,y):\t%.2f %.2f\n", xavg[i], yavg[i]);
+ if( tflag[i] )
+ fprintf(stdout, "\ttol**2:\t\t%.4f (tolerance met)\n", ctolsq);
+ else
+ fprintf(stdout, "\ttol**2:\t\t%.4f\n", ctolsq);
+ fprintf(stdout, "\tregion%d:\t%s\n", i, tbuf);
+ }
+ }
+ /* close the event file */
+ if( !domem ) FunClose(fun);
+ /* close filter for image-style filtering */
+ if( doimfilt || domem ) FilterClose(filter);
+ /* generate a new current region string */
+ *cregion = '\0';
+ /* handle possible first exclude specially */
+ if( regions[0] ){
+ strncat(cregion, regions[0], clen-1);
+ if( cregion[strlen(cregion)-1] != ';' )
+ strncat(cregion, ";", clen-1);
+ }
+ /* handle all other regions */
+ for(i=1; i<=nregion; i++){
+ snprintf(tbuf, SZ_LINE-1, regions[i], xavg[i], yavg[i]);
+ while( (int)(strlen(tbuf) + strlen(cregion)) >= (clen-2) ){
+ clen += SZ_LINE;
+ cregion = xrealloc(cregion, clen);
+ }
+ strncat(cregion, tbuf, clen-1);
+ if( (i != nregion) && (cregion[strlen(cregion)-1] != ';') )
+ strncat(cregion, ";", clen-1);
+ }
+ /* check for no tolerances left to do */
+ if( ntol == 0 ){
+ if( verbose > 1 ) fprintf(stdout, "All tolerances have been met.\n");
+ break;
+ }
+ else{
+ /* save averages for this iteration for next tolerance check */
+ memcpy(pxavg, xavg, sizeof(double)*(MAXREGIONS+1));
+ memcpy(pyavg, yavg, sizeof(double)*(MAXREGIONS+1));
+ }
+ }
+
+ /* final count of events using current region */
+ if( verbose >= 2 ){
+ if( iter > 0 ) fprintf(stdout, "\n");
+ fprintf(stdout, "final: %s\n", cname);
+ }
+ memset(nev, 0, sizeof(int)*(MAXREGIONS+1));
+ if( !domem ){
+ /* append current region to filename root */
+ if( !doimfilt ){
+ /* if the region string is not long, we append it to the file name */
+ if( (strlen(fname2) + strlen(cregion) + strlen(extn)) < (SZ_LINE-1) ){
+ snprintf(cname, clen-1, "%s%s%s", fname2, cregion, extn);
+ }
+ /* otherwise we have to use a temp file */
+ else{
+ /* delete previous temp region file */
+ if( *cregionfile ) unlink(cregionfile);
+ /* get new temp file name */
+ if( !(cfd=mkrtemp(prefix, ".reg", cregionfile, SZ_LINE, 1)) ){
+ gerror(stderr, "can't open temp region file: %s\n", cregionfile);
+ }
+ else{
+ if( !write(cfd, cregion, strlen(cregion)) ){
+ gerror(stderr, "can't write temp file: %s\n", cregionfile);
+ }
+ close(cfd);
+ }
+ snprintf(cname, clen-1, "%s@%s%s", fname2, cregionfile, extn);
+ }
+ }
+ /* but not for image-style filtering -- we filter the region manually */
+ else{
+ strncpy(cname, fname, clen-1);
+ }
+ /* open the data file with the new region as part of the filter */
+ if( !(fun = FunOpen(cname, "r", NULL)) ){
+ gerror(stderr, "can't FunOpen file (or find extension): %s\n", fname);
+ }
+ /* make sure we have a table */
+ FunInfoGet(fun, FUN_TYPE, &type, 0);
+ if( type != FUN_TABLE ){
+ gerror(stderr,
+ "please specify a table (image support not yet implemented)\n");
+ }
+ got = FunColumnSelect(fun, sizeof(EvRec), NULL,
+ "$x", "D", "r", FUN_OFFSET(Ev, x),
+ "$y", "D", "r", FUN_OFFSET(Ev, y),
+ "$region", "J", "r", FUN_OFFSET(Ev, region),
+ NULL);
+ }
+ if( domem || doimfilt ){
+ FunInfoGet(fun, FUN_ENDIAN, &endian, FUN_HEADER, &header,
+ FUN_SECT_X0, &x0, FUN_SECT_X1, &x1,
+ FUN_SECT_Y0, &y0, FUN_SECT_Y1, &y1,
+ FUN_SECT_BLOCK, &block, FUN_BINCOLS, &bincols,
+ FUN_RAWSIZE, &rawsize,
+ 0);
+ if( doimfilt ){
+ /* make sure we will have a valid image section in which to filter */
+ if( (x1-x0 <= 0) || (y1-y0 <= 0) )
+ gerror(stderr, "invalid or zero image dimensions(s) for table\n",
+ bincols ? bincols : "invalid bincols?");
+ /* make up the mode string */
+ snprintf(mode, SZ_LINE-1,
+ "type=events,convert=%s,evsect=\"%d %d %d %d %d\"",
+ (endian == is_bigendian()) ? "false" : "true",
+ x0, x1, y0, y1, block);
+ }
+ else{
+ snprintf(mode, SZ_LINE-1,
+ "type=events,convert=%s",
+ (endian == is_bigendian()) ? "false" : "true");
+ }
+ /* add columns */
+ if( bincols ){
+ strcat(mode, ",");
+ strcat(mode, bincols);
+ }
+ filter = FilterOpen(header, cregion, mode);
+ }
+
+ /* get wcs info */
+ FunInfoGet(fun, FUN_WCS, &wcs, 0);
+ /* and set output */
+ if( (dowcs = (wcs && iswcs(wcs))) ){
+ wcsoutinit(wcs, wcssys);
+ }
+
+ /* filter events and accumulate total events for each region */
+ while( 1 ){
+ if( !domem ){
+ if( !(ebuf = (Ev)FunTableRowGet(fun, NULL, maxrow, NULL, &got)) )
+ break;
+ }
+ /* standard event filtering */
+ if( !domem && !doimfilt ){
+ for(i=0; i<got; i++){
+ ev = ebuf+i;
+ if( ev->region > nregion ){
+ gerror(stderr, "compound regions are not allowed!\n");
+ }
+ nev[ev->region] += 1;
+ }
+ }
+ /* memory and image-style filtering is done manually */
+ else{
+ /* get pointer to raw buffer */
+ FunInfoGet(fun, FUN_RAWBUF, &rawbuf, 0);
+ /* get events which pass the region filter */
+ if( filter && (filter != NOFILTER) &&
+ FilterEvents(filter, rawbuf, rawsize, got, rbuf)){
+ /* loop through events, process those which are in a region */
+ for(i=0; i<got; i++){
+ if( rbuf[i] > 0 ){
+ ev = ebuf + i;
+ nev[rbuf[i]] += 1;
+ }
+ }
+ }
+ }
+ if( ebuf ) xfree(ebuf);
+ if( domem ) break;
+ }
+
+ /* final display */
+ if( verbose > 1) fprintf(stdout, "\n");
+ for(i=1; i<=nregion; i++){
+ /* if we did not iterate, then final x, y is same as initial x, y */
+ if( niter == 0 ){
+ xavg[i] = ox[i];
+ yavg[i] = oy[i];
+ }
+ if( dowcs ){
+ /* convert physical pixels to image */
+ xim = tlp2i(xavg[i], tlmin[0], binsiz[0], tltyp[0]);
+ yim = tlp2i(yavg[i], tlmin[1], binsiz[1], tltyp[1]);
+ pix2wcs(wcs, xim, yim, &ra, &dec);
+ }
+ if( verbose <= 0 ){
+ fprintf(stdout, "%d", nev[i]);
+ if( nev[i] > 0 ){
+ fprintf(stdout, " %.2f %.2f", xavg[i], yavg[i]);
+ if( dowcs ){
+ fprintf(stdout, " %.6f %.6f %s", ra, dec, wcssys);
+ }
+ }
+ else{
+ fprintf(stdout, " 0.00 0.00");
+ if( dowcs ){
+ fprintf(stdout, " 0.0000 0.0000 none");
+ }
+ }
+ fprintf(stdout, "\n");
+ }
+ else{
+ if( i == 1 ){
+ fprintf(stdout, "event_file:\t%s\n", fname);
+ fprintf(stdout, "initial_region:\t%s\n", iregion);
+ fprintf(stdout, "iterations:\t%d\n", niter);
+ fprintf(stdout, "tolerance:\t%.4f\n", tol);
+ }
+ fprintf(stdout, "\n");
+ fprintf(stdout, "events:\t\t%d\n", nev[i]);
+ if( nev[i] > 0 ){
+ fprintf(stdout, "x,y(physical):\t%.2f %.2f\n", xavg[i], yavg[i]);
+ if( dowcs ){
+ fprintf(stdout, "ra,dec(%s):\t%.6f %.6f\n", wcssys, ra, dec);
+ }
+ snprintf(tbuf, SZ_LINE-1, regions[i], xavg[i], yavg[i]);
+ fprintf(stdout, "final_region%d:\t%s\n", i, tbuf);
+ }
+ }
+ if( regions[i] ) xfree(regions[i]);
+ }
+ FunClose(fun);
+ /* close filter for manual filtering */
+ if( doimfilt || domem ) FilterClose(filter);
+ /* handle overflow region specially */
+ if( regions[0] ) xfree(regions[0]);
+ if( cregion ) xfree(cregion);
+ if( cname ) xfree(cname);
+ if( rbuf ) xfree(rbuf);
+ if( *cregionfile ) unlink(cregionfile);
+ /* success */
+ return(0);
+}
generated by cgit v0.12 at 2024-08-29 00:09:46 (GMT)