summaryrefslogtreecommitdiffstats
path: root/funtools/funimage.c
diff options
context:
space:
mode:
Diffstat (limited to 'funtools/funimage.c')
-rw-r--r--funtools/funimage.c581
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);
+}
generated by cgit v0.12 at 2024-08-29 01:12:20 (GMT)