diff options
Diffstat (limited to 'funtools/funmerge.c')
-rw-r--r-- | funtools/funmerge.c | 539 |
1 files changed, 539 insertions, 0 deletions
diff --git a/funtools/funmerge.c b/funtools/funmerge.c new file mode 100644 index 0000000..f141716 --- /dev/null +++ b/funtools/funmerge.c @@ -0,0 +1,539 @@ +/* + * Copyright (c) 1999-2003 Smithsonian Astrophysical Observatory + */ + +#include <funtools.h> +#include <fitsy.h> +#include <wcs.h> +#include <gio.h> +#include <tl.h> +#include <word.h> + +/* mode flags for processing */ +#define FLAG_FILE 1 /* save file name for each event */ +#define FLAG_WCS 2 /* do wcs conversion */ +#define FLAG_SAVE 4 /* do wcs conversion and save old event values */ +#define FLAG_IWCS 8 /* use image wcs */ + +#define MAXROW 8192 +static int maxrow=MAXROW; + +extern char *optarg; +extern int optind; + +typedef struct evstruct{ + double x, y; + double ox, oy; + short nfile; +} *Ev, EvRec; + +typedef struct filestruct{ + struct filestruct *next; + char *fname; + Fun fun; + FITSHead header; + struct WorldCoor *wcs; + int type; + int tltyp[2]; + double tlmin[2]; + double tlmax[2]; + double binsiz[2]; + int bin[2]; + char *bincols[2]; + int nfile; +} *MFile, MFileRec; + +#ifdef ANSI_FUNC +static void InputFileParams(MFile ifile, int flag) +#else +static void InputFileParams(ifile, flag) + MFile ifile; + int flag; +#endif +{ + int i; + int got; + + /* get bin offsets for bin columns */ + FunInfoGet(ifile->fun, FUN_BINOFFS, ifile->bin, 0); + for(i=0; i<=1; i++){ + if( ifile->bin[i] >=0 ){ + FunColumnLookup(ifile->fun, NULL, ifile->bin[i], &(ifile->bincols[i]), + &(ifile->tltyp[i]), NULL, NULL, NULL, NULL); + ifile->tlmin[i] = + FunParamGetd(ifile->fun, "TLMIN", ifile->bin[i]+1, 0.0, &got); + ifile->tlmax[i] = + FunParamGetd(ifile->fun, "TLMAX", ifile->bin[i]+1, 0.0, &got); + ifile->binsiz[i] = + FunParamGetd(ifile->fun, "TDBIN", ifile->bin[i]+1, 1.0, &got); + } + } + /* get WCS info if necessary */ + if( flag & FLAG_WCS ){ + if( flag & FLAG_IWCS ) + FunInfoGet(ifile->fun, FUN_WCS, &(ifile->wcs), 0); + else + FunInfoGet(ifile->fun, FUN_WCS0, &(ifile->wcs), 0); + if( !ifile->wcs ) + gerror(stderr, "no WCS found in %s\n", ifile->fname); + } + /* fitsy header is used for consistency check oncolumns */ + FunInfoGet(ifile->fun, FUN_HEADER, &(ifile->header), 0); +} + +#ifdef ANSI_FUNC +static int OpenInputFile(MFile ifile, int flag) +#else +static int OpenInputFile(ifile, flag) + MFile ifile; + int flag; +#endif +{ + char *mode; + + /* mode depends on file count */ + if( ifile->nfile == 1 ) + mode = "rc"; + else + mode = "r"; + + /* open input file */ + if( !(ifile->fun = FunOpen(ifile->fname, mode, NULL)) ){ + gerror(stderr, "can't FunOpen input file (or find extension): %s\n", + ifile->fname); + return 0; + } + + /* check data type */ + FunInfoGet(ifile->fun, FUN_TYPE, &(ifile->type), 0); + if( (ifile->type != FUN_TABLE) && (ifile->type != FUN_EVENTS) ){ + gerror(stderr, "Sorry, only binary tables can be merged (thus far)\n"); + return 0; + } + + /* get params */ + InputFileParams(ifile, flag); + + /* return the good news */ + return 1; +} + +#ifdef ANSI_FUNC +static MFile NewInputFile( MFile head, char *s, int flag) +#else + static MFile NewInputFile(head, s, flag) + MFile head; + char *s; + int flag; +#endif +{ + MFile thead, cfile, cur; + + /* allocate space */ + if( !(cfile = (MFile)xcalloc(1, sizeof(MFileRec))) ){ + gerror(stderr, "can't allocate input file record for %s\n", s); + return NULL; + } + /* fill in the blanks */ + cfile->fname = xstrdup(s); + /* add to list */ + if( head == NULL ){ + thead = cfile; + cfile->nfile = 1; + /* we open the first input file, which the the base file for wcs */ + if( !OpenInputFile(cfile, flag) ){ + gerror(stderr, "can't open input file: %s\n", cfile->fname); + xfree(cfile); + return NULL; + } + } + else{ + thead = head; + for(cur=head; cur->next!=NULL; cur=cur->next) + ; + cur->next = cfile; + cfile->nfile = cur->nfile+1; + } + /* return head of list */ + return thead; +} + +#ifdef ANSI_FUNC +static void FreeInputFile(MFile file) +#else +static void FreeInputFile(file) + MFile file; +#endif +{ + if( file->fname ) xfree(file->fname); + if( file ) xfree(file); +} + +#ifdef ANSI_FUNC +static void +ConvertEvents(MFile ifile, MFile ofile, MFile file1, Ev ebuf, int nrow, + int flag) +#else +static void +ConvertEvents(ifile, ofile, file1, ebuf, nrow, flag) + MFile ifile; + MFile ofile; + MFile file1; + Ev ebuf; + int nrow; + int flag; +#endif +{ + int i; + int offscl; + double dval1, dval2; + Ev ev=NULL; + + for(i=0; i<nrow; i++){ + if( flag & (FLAG_WCS | FLAG_FILE) ){ + ev = ebuf+i; + } + /* do wcs translation, if necessary */ + if( flag & FLAG_WCS ){ + /* save old coord, might have to write it back out */ + ev->ox = ev->x; + ev->oy = ev->y; + /* wcs convert all but the first file (which is the base file) */ + if( ifile != file1 ){ + if( flag & FLAG_IWCS ){ + /* convert physical pixels to image */ + ev->x = tlp2i(ev->x, + ifile->tlmin[0], ifile->binsiz[0], ifile->tltyp[0]); + ev->y = tlp2i(ev->y, + ifile->tlmin[1], ifile->binsiz[1], ifile->tltyp[1]); + } + /* convert image pixels to ra/dec using wcs */ + pix2wcs(ifile->wcs, ev->x, ev->y, &dval1, &dval2); + /* convert ra/dec to image pixels using base wcs */ + wcs2pix(file1->wcs, dval1, dval2, &(ev->x), &(ev->y), &offscl); + if( flag & FLAG_IWCS ){ + /* convert pixels back to physical in base system */ + ev->x = tli2p(ev->x, + file1->tlmin[0], file1->binsiz[0], file1->tltyp[0]); + ev->y = tli2p(ev->y, + file1->tlmin[1], file1->binsiz[1], file1->tltyp[1]); + } + /* if output coords are not float, we need to round before update */ + switch( ofile->tltyp[0] ){ + case 'D': + case 'E': + break; + default: + if( ev->x >=0 ) ev->x += 0.5; + else ev->x -= 0.5; + break; + } + switch( ofile->tltyp[1] ){ + case 'D': + case 'E': + break; + default: + if( ev->y >=0 ) ev->y += 0.5; + else ev->y -= 0.5; + break; + } + } + } + /* save file number, if necessary */ + if( flag & FLAG_FILE ){ + ev->nfile = ifile->nfile; + } + } +} + +#ifdef ANSI_FUNC +static void +usage (char *fname) +#else +static void usage(fname) + char *fname; +#endif +{ + fprintf(stderr, + "usage: %s [-w|-x] [-f colname] iname1 iname2 ... oname\n", fname); + fprintf(stderr, "optional switches:\n"); + fprintf(stderr, " -f # output a column specifying file from which this event came\n"); + fprintf(stderr, " -w # adjust position values using WCS info\n"); + fprintf(stderr, " -x # adjust position values using WCS info and save old values\n"); + 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 i; + int c; + int iargs; + int nrow; + int flag; + int put; + char tbuf[SZ_LINE]; + char onames[2][SZ_LINE]; + char colname[SZ_LINE]; + char *s; + char *oname=NULL; + char *mode=NULL; + char *filemode=NULL, *savemode=NULL, *wcsmode=NULL; + double tlmin0, tlmin1, tlmax0, tlmax1; + GIO gio; + Ev ev, ebuf; + EvRec tev; + Fun ofun; + MFile ifiles=NULL, ofile=NULL, cfile=NULL, tfile=NULL; + + /* exit on gio errors */ + if( !getenv("GERROR") ) + setgerror(2); + + /* get maxrow,if user-specified */ + if( (s=getenv("FUN_MAXROW")) != NULL ) + maxrow = atoi(s); + + /* we want the args in the same order in which they arrived, and + gnu getopt sometimes changes things without this */ + putenv("POSIXLY_CORRECT=true"); + + /* set colname to something silly. it won't actually be added to the + file unless -c <name> is on the command line, because the mode is + NULL by default */ + strcpy(colname, "???"); + strcpy(onames[0], "???"); + strcpy(onames[1], "???"); + + /* init flag: assume we will use image wcs */ + flag = FLAG_IWCS; + + /* process switch arguments */ + while ((c = getopt(argc, argv, "f:wxIT")) != -1){ + switch(c){ + case 'f': + flag |= FLAG_FILE; + filemode = "w"; + strcpy(colname, optarg); + break; + case 'w': + flag |= FLAG_WCS; + wcsmode = "rw"; + break; + case 'x': + flag |= FLAG_SAVE; + savemode = "w"; + flag |= FLAG_WCS; + wcsmode = "rw"; + break; + case 'I': + flag |= FLAG_IWCS; + break; + case 'T': + flag &= ~FLAG_IWCS; + break; + default: + break; + } + } + + iargs = argc - optind - 1; + /* make sure we have something to do */ + if( iargs < 1 ) usage(argv[0]); + + /* process input files */ + for(i=optind; i<optind+iargs; i++){ + if( argv[i][0] != '@' ){ + /* allocate new ifile (and open the first file, i.e. the base) */ + ifiles = NewInputFile(ifiles, argv[i], flag); + } + /* @ sign means this file contains file names */ + else{ + if( !(gio = gopen(&argv[i][1], "r")) ){ + gerror(stderr, "can't open include file: %s\n", argv[i]); + } + while( ggets(gio, tbuf, SZ_LINE) ){ + if( (*tbuf == '\n') || (*tbuf == '#') ) continue; + nowhite(tbuf, tbuf); + /* allocate new ifile (and open the first file, i.e. the base) */ + ifiles = NewInputFile(ifiles, tbuf, flag); + } + gclose(gio); + } + } + + /* final argument is the output file */ + oname = argv[argc-1]; + /* open the output FITS image, preparing to copy input params */ + if( !(ofun = FunOpen(oname, "w", ifiles->fun)) ) + gerror(stderr, "can't FunOpen output file: %s\n", oname); + /* allocate record for the output file */ + ofile = (MFile)xcalloc(1, sizeof(MFileRec)); + ofile->fun = ofun; + ofile->fname = xstrdup(oname); + InputFileParams(ofile, 0); + + /* get old column names from base file, if necessary */ + if( savemode ){ + for(i=0; i<=1; i++){ + if( ifiles->bincols[i] ){ + snprintf(onames[i], SZ_LINE-1, "OLD_%s", ifiles->bincols[i]); + } + } + } + + /* add merge filename parameters */ + for(cfile=ifiles; cfile; cfile=cfile->next){ + snprintf(tbuf, SZ_LINE-1, "MIF%d", cfile->nfile); + FunParamPuts(ofun, tbuf, 0, cfile->fname, "merge input file", 1); + } + + /* adjust the TLMIN/TLMAX of the base file to accommodate all of the + possible position values from this event file. This effectively + enlargens the size of the output image. */ + /* set temp tlmin variables */ + tlmin0 = ifiles->tlmin[0]; + tlmin1 = ifiles->tlmin[1]; + tlmax0 = ifiles->tlmax[0]; + tlmax1 = ifiles->tlmax[1]; + for(cfile=ifiles->next; cfile; cfile=cfile->next){ + if( !OpenInputFile(cfile, flag) ){ + gerror(stderr, "can't open input file: %s\n", cfile->fname); + } + /* make sure stdin was not specified for one of these input files */ + FunInfoGet(cfile->fun, FUN_GIO, &gio, 0); + if( gio->type & GIO_STREAM ){ + gerror(stderr, "a stream can only be used for the first input file\n"); + } + /* make sure we have the same column structure as the base file */ + if( ifiles->header->table->tfields != cfile->header->table->tfields ){ + gerror(stderr, + "%s and %s have a different number of columns (%d, %d)\n", + ifiles->fname, cfile->fname, + ifiles->header->table->tfields, cfile->header->table->tfields); + } + for(i=0; i<cfile->header->table->tfields; i++){ + if( strcasecmp(ifiles->header->table->col[i].name, + cfile->header->table->col[i].name) ){ + gerror(stderr, + "column names for %s and %s differ in column #%d (%s, %s)\n", + ifiles->fname, cfile->fname, + i+1, + ifiles->header->table->col[i].name, + cfile->header->table->col[i].name); + } + else if( (ifiles->header->table->col[i].type!= + cfile->header->table->col[i].type) || + (ifiles->header->table->col[i].n != + cfile->header->table->col[i].n) ){ + gerror(stderr, + "data types for %s and %s differ in column '%s' (%d%c, %d%c)\n", + ifiles->fname, + cfile->fname, + ifiles->header->table->col[i].name, + ifiles->header->table->col[i].n, + ifiles->header->table->col[i].type, + cfile->header->table->col[i].n, + cfile->header->table->col[i].type); + } + } + ev = &tev; + /* this is the smallest x,y we can have */ + ev->x = cfile->tlmin[0]; + ev->y = cfile->tlmin[1]; + ConvertEvents(cfile, ofile, ifiles, ev, 1, flag); + if( ev->x < tlmin0 ){ + tlmin0 = (int)(ev->x-1); + FunParamPutd(ofun, "TLMIN", cfile->bin[0]+1, tlmin0, 7, + "Min. axis value", 1); + } + if( ev->y < tlmin1 ){ + tlmin1 = (int)(ev->y-1); + FunParamPutd(ofun, "TLMIN", cfile->bin[1]+1, tlmin1, 7, + "Min. axis value", 1); + } + /* this is the largest x,y we can have */ + ev->x = cfile->tlmax[0]; + ev->y = cfile->tlmax[1]; + ConvertEvents(cfile, ofile, ifiles, ev, 1, flag); + if( ev->x > tlmax0 ){ + tlmax0 = (int)(ev->x+1); + FunParamPutd(ofun, "TLMAX", cfile->bin[0]+1, tlmax0, 7, + "Max. axis value", 1); + } + if( ev->y > tlmax1 ){ + tlmax1 = (int)(ev->y+1); + FunParamPutd(ofun, "TLMAX", cfile->bin[1]+1, tlmax1, 7, + "Max. axis value", 1); + } + /* we have to close the file now so as not to run out of descriptors */ + FunClose(cfile->fun); + } + + /* in some cases, we an avoid transfering anything into user space, + by not calling FunColumnSelect() on the input file, but we must + explicitly make a select call to output the raw data */ + if( !filemode && !wcsmode ){ + FunColumnSelect(ofun, 0, "merge=update", NULL); + /* flag that we do not convert data on input or output */ + mode = "convert=false"; + } + + /* now we can process the events in all of the input files */ + for(cfile=ifiles; cfile; cfile=cfile->next){ + /* for all but the base file, re-open input file */ + if( cfile != ifiles ){ + if( !OpenInputFile(cfile, flag) ){ + gerror(stderr, "can't open input file: %s\n", cfile->fname); + } + } + /* select columns for processing */ + if( filemode || wcsmode ){ + FunColumnSelect(cfile->fun, sizeof(EvRec), "merge=update", + "$X", "D", wcsmode, FUN_OFFSET(Ev, x), + "$Y", "D", wcsmode, FUN_OFFSET(Ev, y), + onames[0], "D", savemode, FUN_OFFSET(Ev, ox), + onames[1], "D", savemode, FUN_OFFSET(Ev, oy), + colname, "I", filemode, FUN_OFFSET(Ev, nfile), + NULL); + } + /* make the new extension the reference handle for the output file */ + FunInfoPut(ofile->fun, FUN_IFUN0, &cfile->fun, 0); + /* extract events from this file */ + while( (ebuf = FunTableRowGet(cfile->fun, NULL, maxrow, mode, &nrow)) ){ + /* translate events, if necessary */ + ConvertEvents(cfile, ofile, ifiles, ebuf, nrow, flag); + /* write to output file */ + if( (put=FunTableRowPut(ofile->fun, ebuf, nrow, 0, mode)) != nrow ){ + gerror(stderr, "expected to write %d events; only wrote %d\n", + nrow, put); + } + xfree(ebuf); + } + /* close all files but the first */ + if( cfile != ifiles ) FunClose(cfile->fun); + } + + /* close output before input, so end flush is done automatically */ + FunClose(ofun); + FreeInputFile(ofile); + + /* close base file */ + FunClose(ifiles->fun); + + /* clean up */ + for(cfile=ifiles; cfile;){ + tfile = cfile->next; + FreeInputFile(cfile); + cfile = tfile; + } + return(0); +} |