diff options
Diffstat (limited to 'funtools/funimage.c')
-rw-r--r-- | funtools/funimage.c | 581 |
1 files changed, 581 insertions, 0 deletions
diff --git a/funtools/funimage.c b/funtools/funimage.c new file mode 100644 index 0000000..e763510 --- /dev/null +++ b/funtools/funimage.c @@ -0,0 +1,581 @@ +/* + * Copyright (c) 1999-2005 Smithsonian Astrophysical Observatory + */ + +#include <funtoolsP.h> + +extern char *optarg; +extern int optind; + +#define FORCE_IMAGE -1 + +/* list definitions */ +#define X_DEFAULT "X" +#define Y_DEFAULT "Y" +#define VAL_DEFAULT "VAL" +#define COL1_DEFAULT "COL1" +#define COL2_DEFAULT "COL2" +#define COL3_DEFAULT "COL3" + +/* list event structure */ +typedef struct evstruct{ + double x, y, val; +} *Ev, EvRec; + +/* list structure */ +typedef struct List { + Fun fun; + char *colname[3]; + double tlmin[2]; + double tlmax[2]; + double binsiz[2]; + int dim[2]; + double tscale[2]; + double tzero[2]; + int scaled[2]; +} List; + +#ifdef ANSI_FUNC +static void +doimproj(int dim, int bitpix, void *buf, int i, int j, void *pbuf, int idx) +#else +static void +doimproj(dim, bitpix, buf, i, j, pbuf, idx) + int dim; + int bitpix; + void *buf; + int i; + int j; + void *pbuf; + int idx; +#endif +{ + switch(bitpix/FT_WORDLEN){ + case TY_CHAR: + ((unsigned char *)pbuf)[idx] += ((unsigned char *)buf)[(j*dim)+i]; + break; + case TY_USHORT: + ((unsigned short *)pbuf)[idx] += ((unsigned short *)buf)[(j*dim)+i]; + break; + case TY_SHORT: + ((short *)pbuf)[idx] += ((short *)buf)[(j*dim)+i]; + break; + case TY_INT: + ((int *)pbuf)[idx] += ((int *)buf)[(j*dim)+i]; + break; + case TY_LONG: +#if HAVE_LONG_LONG == 0 + gerror(stderr, + "64-bit data support not built (long long not available)\n"); +#endif + ((longlong *)pbuf)[idx] += ((longlong *)buf)[(j*dim)+i]; + break; + case TY_FLOAT: + ((float *)pbuf)[idx] += ((float *)buf)[(j*dim)+i]; + break; + case TY_DOUBLE: + ((double *)pbuf)[idx] += ((double *)buf)[(j*dim)+i]; + break; + default: + break; + } +} + +#ifdef ANSI_FUNC +static void +usage (char *fname) +#else +static void usage(fname) + char *fname; +#endif +{ + fprintf(stderr, "usage:\n"); + fprintf(stderr, + "%s [-a] iname oname [bitpix=n[,bscale=false]]\n", + fname); + fprintf(stderr, + "%s -l [-a] iname oname xcol:xdims ycol:ydims vcol bitpix=n\n", + fname); + fprintf(stderr, "optional switches:\n"); + fprintf(stderr, " -a # append to existing output file as an image extension\n"); + fprintf(stderr, " -l # input is a list file containing xcol, ycol, value\n"); + fprintf(stderr, " -p [x|y] # project along x or y axis to create a 1D image\n"); + fprintf(stderr, "\n"); + fprintf(stderr, "Specify output image data type using bitpix=n, where n is:\n"); + fprintf(stderr, " 8 # unsigned char\n"); + fprintf(stderr, " 16 # short\n"); + fprintf(stderr, " 32 # int\n"); + fprintf(stderr, " -32 # float\n"); + fprintf(stderr, " -64 # double\n"); + fprintf(stderr, "(Defaults: same as input image, or 32 (int) for tables.)\n"); + fprintf(stderr, "\n"); + fprintf(stderr, "List format requires x, y column information of the form:\n"); + fprintf(stderr, " name:dim # values range from 1 to dim\n"); + fprintf(stderr, " name:min:max # values range from min to max\n"); + fprintf(stderr, " name:min:max:binsiz # dimensions scaled by binsize\n"); + fprintf(stderr, "\n(version: %s)\n", FUN_VERSION); + exit(1); +} + +#ifdef ANSI_FUNC +static List * +ColumnInfo(Fun fun, char *colstr[]) +#else +static List * +ColumnInfo(fun, colstr) + Fun fun; + char *colstr[]; +#endif +{ + int i; + int ip; + char tbuf[SZ_LINE]; + List *list; + + /* allocate list struct */ + if( !(list = (List *)xcalloc(sizeof(List), 1)) ){ + gerror(stderr, "can't allocate list record\n"); + return NULL; + } + + /* process column info */ + newdtable(":"); + for(i=0; i<3; i++){ + ip = 0; + /* column name */ + if( !*colstr[i] || (*colstr[i] == ':') || !word(colstr[i], tbuf, &ip) ){ + switch(i){ + case 0: + if(FunColumnLookup(fun, COL1_DEFAULT, 0, + NULL, NULL, NULL, NULL, NULL, NULL)) + strcpy(tbuf, COL1_DEFAULT); + else{ + strcpy(tbuf, X_DEFAULT); + } + break; + case 1: + if(FunColumnLookup(fun, COL2_DEFAULT, 0, + NULL, NULL, NULL, NULL, NULL, NULL)) + strcpy(tbuf, COL2_DEFAULT); + else{ + strcpy(tbuf, Y_DEFAULT); + } + break; + case 2: + if(FunColumnLookup(fun, COL3_DEFAULT, 0, + NULL, NULL, NULL, NULL, NULL, NULL)) + strcpy(tbuf, COL3_DEFAULT); + else{ + strcpy(tbuf, VAL_DEFAULT); + } + break; + } + } + list->colname[i] = xstrdup(tbuf); + + /* back up to the previous ':' as required by column dimension routine */ + if( ip && (lastdelim() == ':') ) ip--; + /* column dimensions */ + if( colstr[i][ip] != '\0' ){ + switch(i){ + case 0: + case 1: + _FunColumnDims(&colstr[i][ip], 'I', + &list->tlmin[i], &list->tlmax[i], + &list->binsiz[i], &list->dim[i], + &list->tscale[i], &list->tzero[i], &list->scaled[i]); + if( list->tlmin[i] && list->tlmax[i] ) + list->dim[i] = (int)tldim(list->tlmin[i], list->tlmax[i], + list->binsiz[i], 'I'); + else + list->dim[i] = 0; + break; + case 2: + break; + } + } + /* must have dimesion info for x and y */ + else{ + switch(i){ + case 0: + case 1: + gerror(stderr, "dimension specification missing for %s\n", colstr[i]); + return NULL; + case 2: + break; + } + } + } + freedtable(); + + /* set up to read the specified columns */ + FunColumnSelect(fun, sizeof(EvRec), NULL, + list->colname[0], "D", "r", FUN_OFFSET(Ev, x), + list->colname[1], "D", "r", FUN_OFFSET(Ev, y), + list->colname[2], "D", "r", FUN_OFFSET(Ev, val), + NULL); + + /* return list */ + return list; +} + +#ifdef ANSI_FUNC +int +main (int argc, char **argv) +#else +int +main(argc, argv) + int argc; + char **argv; +#endif +{ + int c; + int i, j; + int args; + int dim1, dim2; + int idx; + int dim=0; + int maxidx=0; + int type=0; + int bitpix=0; + int doscale=1; + int doproj=0; + char *s, *t; + char *mode=NULL; + char *tmode=NULL; + char *omode="w"; + char *iname=NULL; + char *oname=NULL; + char *colstr[3]; + char *buf=NULL; + char *pbuf=NULL; + char tbuf[SZ_LINE]; + char newname[SZ_LINE]; + Fun fun, fun2; + + /* list variables */ + int dolist=0; + int ix, iy; + int got; + List *list=NULL; + EvRec evrec; + + /* 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"); + + /* process switch arguments */ + while ((c = getopt(argc, argv, "ailp:")) != -1){ + switch(c){ + case 'a': + omode = "a"; + case 'i': + type = FORCE_IMAGE; + break; + case 'l': + dolist = 1; + break; + case 'p': + doproj = *optarg; + break; + } + } + + /* sanity checks */ + if( doproj ){ + if( dolist ) gerror(stderr, "-p and -l are mutually exclusive\n"); + switch(doproj){ + case 'x': + case 'y': + break; + case 'X': + case 'Y': + doproj = tolower(doproj); + break; + default: + gerror(stderr, "proj specifier must be 'x' or 'y'\n"); + break; + } + } + + /* process arguments */ + args = argc - optind; + if( args < 2 ) usage(argv[0]); + iname = argv[optind+0]; + oname = argv[optind+1]; + if( dolist ){ + if( args != 6 ) usage(argv[0]); + colstr[0] = argv[optind+2]; + colstr[1] = argv[optind+3]; + colstr[2] = argv[optind+4]; + mode = xstrdup(argv[optind+5]); + } + else{ + if( args > 2 ) mode = xstrdup(argv[optind+2]); + } + + /* check for $NEW and rewrite name as /dev/null[TEXT(),...] */ + if( !strncmp(iname, "$NEW[", 5 ) || !strncmp(iname, "$MASK[", 6 ) ){ + /* for masks, we add to the mode */ + if( !strncmp(iname, "$MASK[", 6 ) ){ + if( mode && *mode ){ + snprintf(tbuf, SZ_LINE-1, "%s,mask=all", mode); + xfree(mode); + mode = xstrdup(tbuf); + } + else{ + mode = xstrdup("mask=all,bitpix=8"); + } + } + /* re-write file name as an empty text table */ + memset(newname, 0, SZ_LINE); + strncpy(newname, "/dev/null[TEXT(x:", SZ_LINE-1); + s = strchr(iname, '[')+1; + t = strchr(s, ','); + strncat(newname, s, t-s); + strncat(newname, ",y:", SZ_LINE-strlen(newname)); + if (t != NULL) { + s = ++t; + t = strchr(s, ','); + strncat(newname, s, t-s); + } + strncat(newname, "),", SZ_LINE-strlen(newname)); + if (t != NULL) { + s = ++t; + strncat(newname, s, SZ_LINE-strlen(newname)); + } + iname = newname; + } + + /* open the input FITS file */ + if( !(fun = FunOpen(iname, "rc", NULL)) ){ + gerror(stderr, "can't FunOpen input file (or find extension): %s\n", + iname); + } + + /* open the output FITS image, preparing to copy input params */ + if( !(fun2 = FunOpen(oname, omode, fun)) ){ + gerror(stderr, "can't FunOpen output file: %s\n", oname); + } + + /* list: read x, y, val from list file */ + if( dolist ){ + /* make sure bitpix=n is specified */ + tmode = xstrdup(mode); + if( _FunKeyword(tmode, "bitpix", NULL, tbuf, SZ_LINE) ){ + bitpix = atoi(tbuf); + } + else{ + gerror(stderr, "missing bitpix specification: '%s'\n", mode); + } + if( tmode ) xfree(tmode); + + /* parse column string */ + if( !(list = ColumnInfo(fun, colstr)) ){ + gerror(stderr, "can't get column info\n"); + } + + /* allocate space for the image */ + if( list->dim[1] ) + maxidx = list->dim[0]*list->dim[1]; + else + maxidx = list->dim[0]; + if( !(buf=(char *)xcalloc(maxidx, ft_sizeof(bitpix))) ){ + gerror(stderr, "can't allocate image buffer of size %d\n", maxidx); + } + + /* process events */ + while( FunTableRowGet(fun, &evrec, 1, NULL, &got) && got){ + /* get x, y indices */ + ix = (int)((evrec.x - list->tlmin[0]) / list->binsiz[0]); + if( list->dim[1] ) + iy = (int)((evrec.y - list->tlmin[1]) / list->binsiz[1]); + else + iy = 0; + /* make sure this pixel is within the image */ + idx = iy*list->dim[0]+ix; + if( (idx < 0) || (idx >= maxidx) ){ + gerror(stderr, + "indices (%d,%d) outside of image array (%d,%d) [tlmin problem?]\n", + ix, iy, list->dim[0], list->dim[1]); + } + /* add pixel value */ + switch(bitpix){ + case 8: + ((unsigned char *)buf)[idx] += (unsigned char)evrec.val; + break; + case 16: + ((short *)buf)[idx] += (short)evrec.val; + break; + case 32: + ((int *)buf)[idx] += (int)evrec.val; + break; + case -32: + ((float *)buf)[idx] += (float)evrec.val; + break; + case -64: + ((double *)buf)[idx] += evrec.val; + break; + } + } + /* write the output image, updating the FITS header from the orig file */ + if( !FunImagePut(fun2, buf, list->dim[0], list->dim[1], bitpix, NULL) ){ + gerror(stderr, "could not FunImagePut: %s\n", argv[2]); + } + } + /* no list: bin input data into an image */ + else{ + /* get required information from funtools structure */ + FunInfoGet(fun, + FUN_SECT_BITPIX, &bitpix, + FUN_SECT_DIM1, &dim1, FUN_SECT_DIM2, &dim2, 0); + if( type != FORCE_IMAGE ) FunInfoGet(fun, FUN_TYPE, &type, 0); + /* set dim2 to 1 for a 1D image */ + if( dim2 == 0 ){ + /* oops, can't project a 1D image */ + if( doproj ) + gerror(stderr, "can't project a 1D image\n"); + else + dim2 = 1; + } + /* dimension of projected 1D image */ + if( doproj ){ + switch(doproj ){ + case 'x': + dim = dim1; + break; + case 'y': + dim = dim2; + break; + } + /* allocate buffer for 1D projected image */ + pbuf = (char *)xcalloc(dim, ABS(bitpix)/FT_WORDLEN); + } + + /* process specific type of data */ + switch(type){ + case FUN_IMAGE: + case FUN_ARRAY: + /* if input data is int */ + if( bitpix >= 0 ){ + /* start by assuming we don't want to apply bscale ... */ + doscale = 0; + if( mode && *mode ){ + tmode = xstrdup(mode); + /* but if output data is float, we will apply bscale */ + if( _FunKeyword(tmode, "bitpix", NULL, tbuf, SZ_LINE) ){ + if( (i=atoi(tbuf)) < 0 ){ + /* apple bscale */ + doscale=1; + } + } + if( tmode ) xfree(tmode); + } + } + /* for float data, we apply bscale (won't actually happen) */ + else{ + doscale = 1; + } + if( !doscale ){ + if( mode && *mode ){ + mode = xrealloc(mode, strlen(mode)+SZ_LINE); + strcat(mode, ",bscale=false"); + } + else{ + mode = xstrdup("bscale=false"); + } + } + switch(doproj ){ + case 0: + buf = NULL; + for(j=1; j<=dim2; j++){ + /* extract and bin 1 row of the data section into an image buffer */ + if( !(buf = FunImageRowGet(fun, buf, j, j, mode)) ) + gerror(stderr, "can't FunImageRowGet: %s\n", iname); + /* write output image, updating the FITS header from the orig file */ + if( !FunImageRowPut(fun2, buf, j, j, 0, 0, 0, NULL) ) + gerror(stderr, "can't FunImageRowPut: %s\n", oname); + } + break; + case 'x': + case 'y': + buf = NULL; + for(j=0; j<dim2; j++){ + /* extract and bin 1 row of the data section into an image buffer */ + if( !(buf = FunImageRowGet(fun, buf, j+1, j+1, mode)) ) + gerror(stderr, "can't FunImageRowGet: %s\n", iname); + /* project each pixel of each line */ + for(i=0; i<dim1; i++){ + idx = (doproj == 'x') ? i : j; + doimproj(dim, bitpix, (void *)buf, i, 0, (void *)pbuf, idx); + } + } + /* write output image, updating the FITS header from the orig file */ + if( !FunImagePut(fun2, pbuf, dim, 0, bitpix, NULL) ) + gerror(stderr, "can't FunImagePut: %s\n", oname); + break; + } + break; + + case FUN_TABLE: + case FUN_EVENTS: + case FORCE_IMAGE: + switch(doproj ){ + case 0: + /* extract and bin the data section into an image buffer */ + if( !(buf = FunImageGet(fun, NULL, mode)) ) + gerror(stderr, "can't FunImageGet: %s\n", iname); + /* write output image, updating the FITS header from the orig file */ + if( !FunImagePut(fun2, buf, 0, 0, 0, NULL) ) + gerror(stderr, "can't FunImagePut: %s\n", oname); + break; + case 'x': + case 'y': + /* extract and bin the data section into an image buffer */ + if( !(buf = FunImageGet(fun, NULL, mode)) ) + gerror(stderr, "can't FunImageGet: %s\n", iname); + /* project each pixel of each line */ + for(j=0; j<dim2; j++){ + for(i=0; i<dim1; i++){ + idx = (doproj == 'x') ? i : j; + doimproj(dim, bitpix, (void *)buf, i, j, (void *)pbuf, idx); + } + } + /* write output image, updating the FITS header from the orig file */ + if( !FunImagePut(fun2, pbuf, dim, 0, bitpix, NULL) ) + gerror(stderr, "can't FunImagePut: %s\n", oname); + break; + } + break; + + default: + gerror(stderr, "internal error: unknown type\n"); + break; + } + } + + /* we can explicitly flush remaining extensions that need copying ... */ + FunFlush(fun2, "copy=remaining"); + + /* ... since we already flushed, we can close files in any order */ + FunClose(fun); + FunClose(fun2); + + /* free up space */ + if( dolist ){ + /* free up memory */ + for(i=0; i<3; i++){ + if( list->colname[i] ) xfree(list->colname[i]); + } + xfree(list); + } + if( buf ) xfree(buf); + if( pbuf ) xfree(pbuf); + if( mode ) xfree(mode); + + return(0); +} |